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EXECUTIVE  SUMMARY 

The  broad  objectives  of  this  project  are  to  (a)  develop  a  computational 
methodology  for  implementing  the  criterion  to  predict  stable  tearing  along  a  general  path 
in  three-dimensions  and  (b)  develop  and  validate  a  mixed  mode  fracture  criterion  using 
the  initial  experimental  data  base  and  completed  computational  methodology. 

First,  a  fully  functional  three-dimensional  crack  growth  algorithm  has  been 
developed  that  is  capable  of  predicting  crack  growth  along  a  general,  three-dimensional 
surface.  The  algorithm,  designated  CRACK3D,  utilizes  the  pre-  and  post-processing 
capabilities  of  ANSYS  while  performing  all  of  the  crack-front  calculations  internally.  The 
computer  code,  designated  CRACK3D,  developed  as  part  of  this  project  is  fully  capable 
of  handling  the  kinematics  of  general  crack  growth  and  determining  the  stress  and 
deformation  states  during  crack  growth.  In  addition,  the  code  can  be  interfaced  with 
subroutines  that  implement  a  broad  range  of  fracture  criterion  to  predict  the  instant  and 
direction  of  crack  growth.  Finally,  the  simulation  of  mixed-mode  crack  growth  under 
three-dimensional  (3D)  conditions,  such  as  the  growth  of  surface  cracks,  corner  cracks, 
embedded  cracks,  and  cracks  with  a  curved  crack  surface  and/or  a  curved  crack  front 
has  also  been  fully  automated  within  CRACK3D.  Since  a  major  portion  of  the  work 
performed  to  finalize  CRACK3D  involved  the  development  of  a  remeshing  capability  in 
3D,  Part  II  of  this  final  report  discusses  the  computational  aspects  of  the  simulation 
procedure  and  associated  algorithms  implemented  for  simulating  arbitrary  3D  crack 
growth  under  general  loading  conditions.  In  particular,  this  paper  discusses  the 
strategies  used  for  implementing  automatic  re-meshing  of  regions  around  growing  crack 
fronts  in  a  3D  body,  along  with  verification  examples. 

Second,  Section  III  of  this  report  describes  a  combined  experimental- 
computational  effort  using  aluminum  alloy  2024-T351  that  was  performed  to  gain  an 
understanding  of  slant  fracture  events  and  to  provide  insight  for  establishing  a  general 
fracture  criterion.  The  stable  tearing  fracture  experiments  were  performed  previously 
under  (a)  combined  tension-torsion  (nominal  mixed-mode  I/Ill)  loading  and  (b)  under 
nominally  Mode  I  load  using  Arcan  specimens. 

Two  types  of  finite  element  models  were  considered  for  the  study  of  slant 
fracture:  (a)  combined  tension-torsion  specimens  containing  stationary,  flat  and  slant 
cracks  subject  to  loads  corresponding  to  the  onset  of  crack  growth,  and  (b)  stable 
tearing  crack  growth  with  slanting  in  a  nominal  Mode  I  Arcan  specimen.  Analysis  results 
reveal  that  there  exists  a  strong  correlation  between  the  direction  of  the  maximum 
effective  plastic  strain  ahead  of  a  crack  front  and  the  orientation  direction  of  slant 
fracture.  In  particular,  it  is  observed  that  (a)  at  the  onset  of  crack  growth,  the  angular 
position  of  the  maximum  effective  plastic  strain  around  the  crack  front  serves  as  a  good 
indicator  for  the  slant  fracture  surface  orientation  during  subsequent  crack  growth;  and 
(b)  during  stable  tearing  crack  growth  with  a  flat-to-slant  transition,  the  crack  growth 
path  on  each  section  plane  through  the  thickness  of  a  specimen  coincides  with  the 
angular  position  of  the  maximum  effective  plastic  strain  around  the  crack  front.  The 
results  of  this  study  suggest  that  the  effective  plastic  strain  may  be  used  as  a  fracture 
parameter  for  predicting  the  fracture  surface  path  of  stable  crack  growth  with  slant 
fracture. 
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I.  INTRODUCTION 

As  noted  by  the  National  Materials  Advisory  Board,  a  critical  limitation  in  the 
state-of-the-art  for  crack-growth  simulation  technology  is  the  lack  of  an  experimentally 
verifiable,  general  stable  tearing  criterion  that  can  be  applied  to  the  whole  spectrum  of 
possible  crack-tip  loading  conditions.  This  would  include  plane  stress,  plane  strain  and 
general,  three-dimensional  conditions. 

This  program  was  directed  towards  (a)  developing  3D  simulation  algorithms  that 
can  be  used  to  perform  node-by-node  stable  tearing  analyses  on  aerospace  alloys  and 
(b)  understanding  and  developing  a  reliable  fracture  criterion  for  prediction  of  stable 
tearing  in  3D  structures. 

II.  Computational  Aspects  and  Code  Development  of  Three-dimensional 
Crack  Growth  Simulations 

Abstract 

An  important  task  in  mixed-mode  fracture  analysis  and  prediction  is  the 
simulation  of  crack  growth  under  mixed-mode  conditions.  To  complete  such  a  task,  one 
must  have  (a)  a  computer  code  capable  of  handling  the  kinematics  of  general  crack 
growth  and  determining  the  stress  and  deformation  states  during  crack  growth,  and  (b) 
a  fracture  criterion  that  can  properly  predict  the  instant  and  direction  of  crack  growth.  A 
current  challenge  is  the  simulation  of  mixed-mode  crack  growth  under  three- 
dimensional  (3D)  conditions,  such  as  the  growth  of  surface  cracks,  corner  cracks, 
embedded  cracks,  and  cracks  with  a  curved  crack  surface  and/or  a  curved  crack  front. 

This  section  addresses  item  (a)  in  the  above  discussion  and  describes  the 
computational  aspects  of  a  simulation  procedure  and  associated  algorithms  for 
simulating  arbitrary  3D  crack  growth  under  general  loading  conditions,  which  have  been 
developed  and  successfully  implemented  by  the  authors  in  a  custom,  finite  element 
based,  crack  growth  analysis  and  simulation  code —  CRACK3D.  In  particular,  this  paper 
will  present  strategies  for  automatic  re-meshing  of  regions  around  growing  crack  fronts 
in  a  3D  body,  and  will  discuss  verification  examples. 


11.1.  Introduction 

Residual  strength  and  fatigue  lifetime  prediction  for  damaged  structures  and 
components  has  been  attracting  more  and  more  attention  of  engineers  and  scientists  [1- 
5].  In  many  cases  the  residual  strength  and  lifetime  of  such  kinds  of  structures  and 
components  is  depending  on  the  behavior  of  flaws  in  them  [6,  7].  Thereby,  the  accurate 
prediction  of  three-dimensional  crack  propagation  in  ductile  materials  plays  an  important 
role  in  the  evaluation  of  residual  strength  and  lifetime  of  structures  and  components  with 
cracks  [8-10]. 

Due  to  the  strict  requirement  on  boundary  conditions  and  the  geometry  of  objects, 
currently  the  analytical  solution  is  nearly  impossible  to  be  obtained  in  many  problems  of 
cracks,  especially  for  the  process  of  ductile  crack  propagation  in  three-dimensional 
structures.  Thus,  different  attempts  have  been  made  to  develop  efficient  and  robust 
numerical  simulation  tools  to  predict  the  response  of  components  during  the  process  of 
crack  propagation  [11-15]. 
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To  develop  a  general  three-dimensional  crack  growth  simulation  tool,  there  are 
several  key  issues.  Firstly,  a  three-dimensional  fracture  criterion  is  required  to  be  used 
to  determine  the  onset  and  the  direction  of  crack  growth.  A  lot  of  experimental, 
analytical  and  numerical  studies  have  been  done  on  this  topic  [16-20]  and  the  results 
show  that,  for  different  kinds  of  materials,  such  as  brittle  materials  and  ductile  materials, 
the  behavior  of  failure  has  more  or  less  difference,  especially  for  mixed-mode  fracture. 
For  aluminum  alloy  2024-T3,  experimental  results  show  that  there  is  a  sharp  transition 
of  crack  growth  behavior  from  predominantly  Mode  I  type  to  Mode  II  type  during  the 
variation  of  the  amount  of  mixed  mode  loading  [18].  Clearly,  using  Stress  Intensity 
Factor  (SIF)  [8,  12]  to  predict  the  direction  of  crack  growth  in  ductile  materials  will  result 
in  inconsistent  results  with  the  experimental  results  [18]. 

Secondly,  a  robust  remeshing  algorithm  is  requisite  to  simulate  the  crack  growth  by 
updating  the  finite  element  mesh.  As  one  knows,  the  geometric  topology  of  structures 
can  be  any  type  during  the  evolution  of  crack  surfaces.  To  make  the  simulation  of  crack 
growth  run  smoothly  and  automatically,  the  remeshing  algorithm  should  have  the  ability 
to  handle  any  kind  of  cracks,  such  as  surface  cracks,  through-thickness  cracks  and 
embedded  cracks  [21].  It  is  usually  required  in  fracture  mechanics  that  the  element  size 
in  the  region  near  crack  fronts  is  quite  small  comparing  to  the  element  size  in  the  region 
far  from  the  crack  fronts.  However,  the  crack  front  may  go  anywhere  and  the  dense 
mesh  in  the  previous  position  of  the  crack  front  contributes  little  to  the  quantities  of  fields 
near  the  current  crack  front.  If  the  remeshing  scheme  can  not  loose  the  mesh  in  the 
previous  position  of  the  crack  front,  the  total  number  of  element  will  increase 
significantly  as  the  movement  of  crack  front  [5],  which  causes  the  size  of  finite  element 
model  and  subsequent  CUP  time-consuming  increasing  quickly  during  the  process  of 
crack  growth. 

In  this  paper,  the  authors  focus  on  the  development  of  new  software,  CRACK3D,  for 
three-dimensional  crack  growth  simulation.  It  addresses  the  computational  aspects  of 
the  simulation  procedure  and  associated  algorithms  used  for  simulating  arbitrary  three- 
dimensional  crack  propagation  under  general  loading  conditions  which  have  been 
developed  and  successfully  implemented  into  a  custom,  finite  element  based,  crack 
growth  simulation  code  —  CRACK3D.  The  paper  is  organized  as  follows.  Section  2 
introduces  the  hierarchy  of  the  simulation  tool,  CRACK3D.  Section  3  describes  the 
schemes  of  modification  of  geometry  and  topology  after  each  step  of  crack  growth. 
Section  4  discusses  the  issues  of  cutting  a  local  region  from  the  global  model  for 
remeshing.  Section  5  gives  the  approaches  of  local  remeshing  and  mesh  connection 
between  the  local  regions  and  the  remnant  of  global  model.  Section  6  gives  the 
approach  used  in  CRACK3D  for  data  transformation  from  the  previous  mesh  to  the 
current  mesh.  Section  7  employs  two  selected  examples  to  demonstrate  the  capabilities 
of  CRACK3D  in  simulating  three-dimensional  crack  growth.  Section  8  presents  the 
conclusions  of  the  work. 

11.2.  Hierarchy  of  the  software  CRACK3D 

Crack  growth  simulation  is  a  combined  incremental  process  including  increasing  or 
decreasing  loads  and  incrementally  variation  of  the  size  of  crack  surface,  which  is  more 
complicated  than  many  other  numerical  applications  in  computational  mechanics  in  the 
aspects  of  the  structure  of  program  and  the  management  of  database.  To  take  the 
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advantage  of  the  existing  pre-  and  post-processing  capabilities  of  general  commercial 
software,  such  as  ANSYS,  current  version  of  CRACK3D  was  developed  to  focus  on  the 
three-dimensional  crack  growth  simulation.  During  the  whole  process  of  crack  growth 
simulation,  CRACK3D  has  no  data-exchange  with  third-party  software.  The  third-party 
software  has  nothing  to  do  with  the  analysis  of  mechanics  except  providing  initial  mesh 
and  visualization  of  results.  Fig.  1  shows  the  relation  of  CRACK3D  with  commercial 
software  ANSYS. 

CRACK3D  is  a  package  of  three-dimensional  crack  propagation  simulation,  which 
was  developed  modularly.  Generally  speaking,  CRACK3D  consists  of  four  modules, 
stress  analysis,  determination  of  new  crack  front,  modification  of  geometry  and  topology 
and  local  remeshing.  These  modules  are  locally  independent  and  invoked  successively 
during  the  simulation  of  crack  growth,  as  shown  in  Fig.  2. 

It  has  been  widely  accepted  that  the  macroscopic  fracture  criteria  are  dependent  on 
the  behavior  and  properties  of  materials,  such  as  brittle  and  ductile  materials.  For  the 
consideration  of  the  implementation  of  different  fracture  criteria  into  CRACK3D  in  the 
future,  a  user  subroutine  interface  is  provided  for  users  to  implement  their  own  fracture 
criteria  into  CRACK3D.  From  Fig.  2,  one  can  find  that  fracture  criteria  will  work  on  the 
base  of  the  stress  and  strain  results  provided  by  the  module  “STRESS  ANALYSIS”;  it 
will  provide  information  necessary  for  determining  the  location  of  new  crack  fronts  and 
the  shape  of  new  crack  surfaces  in  the  module  “NEW  CRACK  FRONT”.  Another  use  of 
the  user  subroutine  interface  of  fracture  criteria  is  that,  users  can  take  the  advantage  in 
studying  or  verifying  new  fracture  criteria  for  different  materials. 

The  module  “GEOMETRY  &  TOPOLOGY  MODIFICATION”  utilizes  the  information 
of  new  crack  front  and  new  crack  surfaces  to  modify  the  configuration  of  the  structure 
under  consideration,  including  its  geometry  and  topology,  and  creates  a  reasonable 
geometric  model  with  the  combination  of  new  crack  surfaces.  The  module  “LOCAL 
REMESHING”  will  generate  a  new  mesh  around  the  new  crack  front  and  convert  the 
numerical  results  from  previous  mesh  to  current  mesh  so  that  the  model  “STRESS 
ANALYSIS”  can  accept  the  data  of  the  new  mesh  and  automatically  continue  the 
simulation. 

It  need  mention  that,  even  though  the  remeshing  tool  is  available  in  current  version 
of  CRACK3D,  the  initial  mesh  is  required  for  user  to  input,  which  can  also  be  generated 
by  a  third-party  software.  Once  the  initial  mesh  data  and  controlling  parameters  are 
provided,  CRACK3D  can  automatically  perform  fracture  analysis,  modification  of 
geometry  and  topology,  remeshing  and  data  transformation  in  sequence  without  user’s 
interference  during  the  whole  process  of  crack  growth  simulation. 

11.3.  Modification  of  geometry  and  topology 

To  simulate  the  crack  growth  in  three-dimensional  structures,  a  fracture  criterion  is 
requisite  to  determine  the  onset  and  direction  of  crack  propagation  at  each  node  along 
the  current  crack  front.  The  distance  of  crack  extension  at  different  nodes  may  be 
assumed  to  be  proportional  to  the  magnitude  of  driving  force  at  the  corresponding  node. 
By  utilizing  the  direction  and  the  length  of  crack  growth  at  each  node  along  crack  front, 
one  can  straightforwardly  obtain  the  coordinates  of  the  location  of  new  crack  front. 

It  is  apparent  that  the  new  crack  front  determined  by  the  direction  and  length  of 
crack  extension  is  a  virtual  crack  front,  which  needs  some  modification  so  that  it 
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becomes  the  real  crack  front.  Fig.  3  illustrates  some  kinds  of  the  modifications  to  the 
virtual  crack  front  for  a  plate  with  a  through-thickness  crack.  In  most  cases,  for 
instance,  some  parts  of  the  virtual  crack  front  may  go  outside  the  material,  especially  for 
a  curvilinear  crack  front.  Firstly,  the  module  “GEOMETRY  &  TOPOLOGY 
MODIFICATION”  determines  the  real  crack  front  by  making  use  of  the  virtual  crack  front 
and  geometric  information  of  the  boundary  surfaces  of  the  structure.  Secondly,  it 
performs  the  geometric  Boolean  calculation  by  utilizing  the  virtual  new  crack  surfaces  to 
determine  the  new  topology  of  the  structure,  and  subsequent  new  boundary  lines,  new 
boundary  surfaces.  Thirdly,  some  of  the  existing  boundary  surfaces  and  boundary  lines 
also  need  to  be  updated  by  the  module  to  reflect  the  change  of  geometry  induced  by  the 
crack  propagation. 

In  order  to  generate  and  provide  a  reasonable  geometric  model  for  the  purpose  of 
remeshing,  some  criteria  are  used  during  the  modification.  For  example,  for  a  through¬ 
thickness  crack  as  shown  in  Fig.  3,  if  the  predicted  location  of  any  end  of  the  crack  front 
is  too  close  to  any  one  of  the  surfaces  of  the  structure  (see  Fig.  3d)  the  ill-shaped 
elements  will  be  generated  in  the  vicinity  of  the  end  of  crack  front  during  remeshing.  To 
solve  this  problem,  a  straightforward  method  is  to  move  the  predicted  location  of  the 
end  of  crack  front  onto  the  closest  surface  to  it.  In  CRACK3D,  a  tolerance  less  than  or 
equal  to  one  fourth  of  the  minimum  element  size  is  used  to  determine  if  the  movement 
is  necessary.  Generally  speaking,  to  avoid  the  occurrence  of  ill-shaped  elements,  the 
modification  to  the  virtual  crack  front  should  result  in  a  real  crack  front  which  satisfies 

1.  There  is  no  node  on  the  crack  front  whose  distance  to  its  closest 
boundary  surface  is  less  than  the  specified  tolerance  and  greater 
than  zero. 

2.  The  length  of  any  one  of  crack  fronts  is  greater  than  the  minimum 
element  size. 

3.  For  surface  cracks  and  through-thickness  cracks,  there  is  no  end  of 
the  crack  fronts  whose  distance  to  its  closest  boundary  line  is  less 
than  the  specified  tolerance  and  greater  than  zero. 

4.  For  each  crack  front,  the  number  of  nodes  on  the  crack  front  whose 
distance  to  its  corresponding  closest  boundary  surface  is  equal  to 
zero  can  only  be  equal  to  zero  or  two. 

Once  the  locations  of  crack  fronts  have  been  determined,  the  module  “GEOMETRY 
&  TOPOLOGY  MODIFICATION”  is  to  perform  the  geometric  Boolean  calculation  by 
utilizing  the  virtual  new  crack  surfaces  to  determine  the  new  topology  of  the  structure, 
and  subsequent  new  boundary  lines,  new  boundary  surfaces.  To  reduce  the  possibility 
of  generating  ill-shaped  elements  during  the  process  of  remeshing,  the  resulting 
geometric  model  should  also  satisfy 

1.  There  is  no  boundary  line  whose  length  is  less  than  the  minimum 
element  size. 

2.  There  is  no  boundary  surface  whose  area  is  equal  to  zero. 

3.  The  area  of  any  separated  cross-section  cut  by  any  oriented  plane  is 
greater  than  a  specified  tolerance. 

4.  Any  segment  of  boundary  lines  must  be  shared  and  can  only  be 
shared  by  two  boundary  surfaces. 
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11.4.  Determination  of  local  region 

To  take  into  account  the  effect  of  geometric  change  induced  by  crack  propagation 
on  the  finite  element  model,  CRACK3D  has  the  capacity  to  specify  a  certain  region  by 
making  use  of  the  user-inputted  parameter  controlling  the  size  of  a  local  region  and 
remesh  it.  A  straightforward  way  is  to  investigate  each  element  to  see  if  the  minimum 
distance  of  its  centroid  to  any  node  on  the  new  crack  front  is  less  than  the  given 
distance  specified  by  users.  A  local  region  will  be  formed  by  collecting  all  the  elements 
lying  within  the  distance.  The  specific  approach  is  given  in  the  following 

1.  Find  out  the  geometric  topology  of  region  Q  on  the  base  of  the 
current  mesh.  Suppose  Sn  is  the  set  of  all  boundary  surfaces  of  Q , 

Ln  is  the  set  of  all  boundary  lines  of  all  boundary  surfaces  of  Q 

sn={s„,s,  sn} 

Ln  ~  | fj)  A  e 

2.  Determine  Q,  (k-1,  2,  ...  n),  which  is  a  local  region  around  the  k?h 
crack  front,  the  size  of  Clk  is  controlled  by  the  parameter  R,  where  n 
is  the  total  number  of  crack  fronts  in  the  structure. 

3.  Determine  AR,  which  is  the  remnant  of  region  Q  subtracted  with  Qk 

(k=1,  2,  ...  n),  thus 

Q=A*U  fi,U  ^2U  ...U  Fln 

4.  Find  out  the  geometric  topology  of  region  AR  based  on  the  current 
mesh  and  the  sets  of  Ln  and  Sn.  Suppose  SA  is  the  set  of  all 
boundary  surfaces  of  AR,  LA  is  the  set  of  all  boundary  lines  of  A  R 

Sa={5„5,eA„} 

A\  —  Ifi’  h  e  ^  r  } 

5.  Investigate  the  relationship  of  all  the  local  regions  Cik  (k=1,  2,  ...  n). 

Check  if  one  of  them  intersects  with  another.  If  yes,  then  combine 
those  local  regions  intersected  mutually  to  form  a  larger  new  local 
region.  Suppose  there  are  m  local  regions  9 \k  ( k-1 ,  2,  ...  m ) 

constructed  after  the  combination  among  the  set  of  Q(  (i=1,  2,  ...n), 
accordingly 

n=AR[j  9i,  u  9?2  U  •••  U  9im 

In  most  cases,  the  local  region  formed  in  the  above  approach  is  acceptable  and  can 
be  put  back  into  global  model  after  remeshing.  We  also  experienced  difficulties  in  some 
cases  that  the  new  mesh  of  the  local  region  formed  in  this  way  after  remeshing  did  not 
match  the  previous  mesh  outside  the  local  region  on  the  interface.  By  studying  the 
findings  we  found  that  the  modification  to  the  local  region  determined  in  the  above 
method  became  requisite  so  that  the  local  region  after  modification  should  meet  the 
following  conditions 
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5.  There  is  no  segment  of  any  global  boundary  lines  in  the  local  region 
S.H/C  (k=1,  2,  ...  m)  which  is  shared  by  the  boundary  line  in  the 

remnant  AR 

6.  For  any  one  of  local  regions  s.R/f  (k=1,  2,  ...m),  there  is  no  segment 
of  the  local  boundary  lines  which  is  shared  by  more  than  two  local 
boundary  surfaces  in  the  local  region  9^  (k=1,  2,  ...  m). 

It  needs  to  be  mentioned  that,  after  each  modification,  the  evaluation  of  the 
relationship  among  all  the  local  regions  and  the  geometric  features  of  each  local  region 
should  be  performed  again  until  no  modification  is  needed  and  all  the  local  regions 
satisfy  the  above  conditions. 

Clearly,  if  the  user-inputted  parameter  controlling  the  size  of  a  local  region  is  large 
enough  to  let  all  the  elements  lie  inside  the  local  region,  then  the  local  region  is  actually 
equal  to  the  whole  global  region,  and  all  the  above  requirements  on  the  local  region  will 
be  met  automatically.  It  will  simplify  the  determination  of  local  region,  but  the  time 
consumed  during  the  process  of  remeshing  and  data  transformation  from  previous 
mesh  to  new  mesh  will  increase  significantly. 

11.5.  Remeshing  on  the  specified  local  region 

Once  the  local  regions  have  been  determined,  the  module  “LOCAL  REMESHING” 
will  be  invoked  to  re-mesh  the  local  region  one  by  one  and  perform  the  data 
transformation  from  previous  mesh  to  the  new  generated  mesh. 

As  one  can  understands,  any  kind  of  geometry  shapes  and  topologies  has  the 
possibility  of  occurrence  during  the  process  of  crack  propagation  that  makes  remeshing 
around  crack  fronts  more  difficult  than  many  other  kinds  of  remeshing  in  a  region 
without  cracks,  and  the  robustness  of  the  remeshing  tool  usually  turn  to  decrease.  To 
develop  a  robust  remeshing  tool  for  arbitrary  three-dimensional  structures  with  cracks, 
two  features  of  the  regions  with  cracks  require  the  most  attention.  One  feature  arises 
due  to  the  geometric  coincidence  of  two  crack  surfaces.  When  this  occurs,  two  nodes 
on  separate  fracture  surfaces  may  occupy  the  same  spatial  position.  Since  there  are 
two  candidate  nodes  at  the  same  spatial  position,  simple  geometric  arguments  are  not 
sufficient  to  identify  the  appropriate  node  for  inclusion  in  a  new  element’s  definition 
during  the  meshing  (re-meshing)  process.  The  second  feature  is  introduced  when  there 
is  inter-penetration  of  two  crack  surfaces,  a  common  occurrence  during  the  meshing 
(re-meshing)  of  curved  crack  surfaces  that  results  in  the  intersection  of  volume 
elements  in  the  region  adjacent  to  curved  crack  surfaces. 

There  are  three  kinds  of  basic  cracks  in  three-dimensional  structures,  surface 
cracks,  through-thickness  cracks  and  embedded  cracks,  which  result  in  different  kinds 
of  geometric  topologies.  The  module  “LOCAL  REMESHING”  in  CRACK3D  was 
developed  based  on  the  mesh  generation  algorithm  given  in  [21]  and  took  into 
consideration  all  the  features  of  the  three  kinds  of  cracks.  Figs  4-6  show  the  resulting 
meshes  of  domains  containing  different  kinds  of  cracks.  From  the  resulting  meshes,  it 
can  be  found  that  the  element  quality  around  any  kind  of  crack  fronts  is  acceptable. 

As  mentioned  in  section  4,  if  a  local  region  is  to  be  remeshed  then  it  will  become 
more  complicate  than  remeshing  the  global  region  since  there  exists  a  necessary 
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requirement  during  the  local  remeshing.  It  is  clear  that,  since  only  parts  of  the  global 
region  undergoes  remeshing,  the  mesh  consistency  on  the  interface  between  the  local 
regions  and  the  remnant  of  the  global  region  obtained  by  subtracting  the  local  regions 
should  be  satisfied  so  that  the  local  regions  after  remeshing  can  be  put  back  and 
connected  to  the  previous  mesh  around  the  local  regions  without  any  gap  and/or 
penetration.  This  process  can  be  illustrated  in  Fig.  7. 

To  assure  the  consistency  of  mesh  on  the  interface  between  local  regions  and  the 
remnant  of  the  global  region,  it  needs  to  identify  the  local  boundary  lines  and  local 
boundary  surfaces  which  need  to  be  remeshed,  and  those  which  should  not  be 
remeshed  so  that  the  profile  of  the  local  mesh  on  the  interface  keep  consistent  with  that 
of  the  remnant  mesh.  To  achieve  this  goal,  the  following  approach  is  suggested 

1.  For  a  given  global  region  Q,  denote  Sn  the  set  of  all  boundary 
surfaces  of  Cl,  Ln  the  set  of  all  boundary  lines  of  all  boundary 
surfaces  of  Q 

Sn  =  {s„S,e£l} 

—  {/,  ;  /;  e 

2.  For  any  one  of  local  regions,  9?/(  (k-1,  2,  ...  m ),  denote  the  set 
of  all  local  boundary  surfaces  of  the  local  region  9^ ,  L)h  the  set  of 
all  local  boundary  lines  of  the  local  region  9J* , 

S,,=|s„Sts*,l  (k=1,2,...m) 

(k=1,  2,  ...  m) 

3.  Denote  A  R  the  remnant  of  the  global  region  Cl  subtracted  with  9t* 

(k-1,  2,  ...  m),  which  satisfies 

Q=Ar  U  %  U  9t2  U  ...U  SR* 

4.  Evaluate  the  geometric  topology  of  the  remnant  region,  AR. 
Suppose  SA  is  the  set  of  all  boundary  surfaces  of  A R,  LA  is  the  set 
of  all  boundary  lines  of  all  boundary  surfaces  of  A  R 

SA={5„5,eA„} 

A\  —  {A  h  e  } 

5.  Evaluate  the  relationship  between  SA  and  S.}ij  (k=1,  2,  ...  m ),  and 
the  relationship  between  LA  and  Lm  (k=1,  2,  ...  m).  Denote  Sc  the 
common  part  of  SA  and  S.}h  (k=1,  2,  ...  m),  Lc  the  common  part  of 
La  and  Z3i(  (k=1,  2,  ...  m) 

3c~^a  n  n  s.ju  n  •••(!  s.Kt 
lc~la  n  n  l,j{i  n  •••/!  Ait, 

6.  For  any  one  of  local  regions,  9?^  (k=1,  2,  ...  m),  decompose  and 

into  two  parts,  S'Rt  and  S\t ,  l}R  and  L2R  ,  respectively,  which 
yield  to  the  following  relations 
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S«,=Si,USj,,  Sj,  «SC,  S-sS,  (k=1,  2,  ...  m) 
knt  —  LR)  U  LRt ,  LRk  £  Lc ,  Lr  e  Lc ,  (k—1,  2,  ...  m) 

7.  For  the  local  region  9^  (k=1,  2,  ...  m),  remesh  the  boundary  lines  in 
L]Ri  and  the  boundary  surfaces  in  SR  ,  while  the  mesh  on  the 
boundary  lines  in  l?Rt  and  on  the  boundary  surfaces  in  keep  the 

original.  After  all  the  mesh  on  the  boundary  lines  and  the  boundary 
surfaces  of  the  local  region  9 \k  is  determined,  the  volume  mesh  can 

be  generated  in  the  local  region,  9?/( . 

After  remeshing  in  each  local  region,  the  techniques  of  mesh  optimization  are 
employed  in  the  module  “LOCAL  REMESHING”  to  optimize  the  distribution  and  quality 
of  elements  around  the  crack  fronts.  Since  the  gradient  of  element  size  in  the  region 
near  crack  fronts  is  severe,  the  resulting  meshes  after  mesh  optimization  also  show  that 
the  transition  of  element  size  gets  significant  improvement  in  such  regions. 

11.6.  Data  mapping  from  previous  mesh  to  current  mesh 

For  the  incremental  elastic-plastic  stress  analysis,  the  numerical  calculation  in  time 
step  fw  is  based  on  the  numerical  results  in  time  step  t^-i-  To  simulate  the  crack 
propagation  automatically  after  remeshing  without  interruption  due  to  the  evolution  of 
finite  element  model  induced  by  crack  growth,  the  previous  results,  such  as 
displacement,  strain  and  stress,  need  to  be  mapped  from  the  previous  mesh  to  the 
current  mesh  so  that  the  subsequent  computation  based  on  the  current  mesh  has  valid 
and  continuous  results. 

Currently  there  are  two  types  of  elements  used  in  CRACK3D  for  the  three- 
dimensional  simulation  of  crack  propagation  with  local  remeshing,  4-noded  tetrahedral 
element  and  10-noded  tetrahedral  element.  Hereinafter  we  will  employ  the  4-noded 
tetrahedral  elements  to  demonstrate  the  approach  used  in  CRACK3D  for  the  data 
transformation  between  two  adjacent  mesh  frames. 


11.7.  Locate  a  new  node  in  the  previous  mesh  frame 

Suppose  ( xo ,  yo,  zo)  is  the  coordinates  of  a  given  point  P  in  the  current  mesh  frame, 
which  may  be  either  a  new  nodal  point  or  a  Gaussian  integration  point  of  a  new  element 
in  the  current  finite  element  mesh,  we  need  search  an  element  first  in  the  previous 
mesh  which  contains  the  given  point  P. 

For  an  arbitrary  element,  say  EP,  in  the  previous  mesh,  without  loss  of  generality, 
suppose  Ni,  N2,  N3  and  N4  are  the  four  nodes  of  the  element  EP,  as  shown  in  Fig.  8. 
One  can  calculate  the  volumes  of  the  four  tetrahedrons  PN2N3N4 ,  PN3NlN4 ,  PNlN2N4 

and  PNlN3N2  by  evaluating  the  following  determinants 


T2  "To 
Ts  "To 
T4  "To 
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(la) 


(1b) 


*3“*0  ^3-^0  Z3~Z0 

V2=-xx-x0  yx-y0  zx-z0 

6 

^4-^0  74  "Jo  Z4~Z0 

7i-7o  Z1"Z0 

J2-Jo  Z2  Z0  (1C) 

J4-Jo  Z4  Z0 

*i-*o  Ji-Jo  Z1-Z0 

V4=-x3-x0  y3-y0  Z3-Z0  (Id) 

6 

x2-x0  y2-y  o  z2-z0 

where  (x,,  y,,  zi)  ( i=1 ,  2,  3,  4)  are  the  coordinates  of  nodes  Ni,  N2,  N3  and  N4.  By  utilizing 
the  magnitudes  of  Vi,  V2,  V3,  and  V4,  one  can  draw  a  conclusion  whether  or  not  the 
given  point  P  is  inside  the  element  EP.  If  none  of  the  volumes  Vi,  V2,  V3,  and  V4  is  less 
than  zero,  then  the  element  EP  contains  the  given  point  P,  otherwise  it  doesn’t  contain 
the  given  point  P. 

To  speed  up  the  search  process  in  CRACK3D,  only  the  elements  of  the  previous 
mesh  inside  the  local  region  are  chosen  to  be  the  candidates. 

11.8.  Evaluate  the  field  quantities  for  a  specified  point 

Once  the  element  in  the  previous  mesh  containing  the  point  in  the  current  mesh  has 
been  found,  a  straightforward  approach  is  used  to  evaluate  all  the  field  quantities  such 
as  displacement,  strain  and  stress,  at  the  new  point  in  the  current  mesh. 

For  a  given  new  point  P  in  the  current  mesh,  suppose  the  element  EP  in  the 
previous  mesh  contain  the  given  point  P.  Then  one  can  evaluate  the  field  quantities  at 
the  point  P  by  employing  the  field  quantities  of  element  EP  in  the  frame  of  previous 
mesh.  For  any  specific  field  quantity  at  point  P,  FP,  it  can  be  interpolated  by  making  use 
of  the  shape  function  of  the  4-noded  tetrahedral  element  with  respect  to  the  location  of 
point  P. 

The  shape  function  of  the  4-noded  tetrahedral  element  with  respect  to  the 
coordinates  of  point  P  can  be  expressed  as 

®1=VJVE,  ®2=v2/vE,  o3=v3/ve,  ®4  =  f4/ve, 

where  VE  =  VX  +V2  +V3  +V4.  Accordingly,  the  field  quantity  at  point  P,  FP,  can  be 
interpolated  in  the  following 

4 

F!  = 

1=1 

where  Fkp  is  the  quantity  at  the  point  P  in  mesh  k,  Fk~l  is  the  quantity  at  node  /  of 
element  EP  in  mesh  k-1 . 
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11.9.  Geometry  and  Topology  Changes  due  to  Crack  Evolution 

The  module  “GEOMETRY  &  TOPOLOGY  MODIFICATION”  in  CRACK3D  can 
handle  geometric  and  topological  changes  due  to  several  common  types  of  crack 
evolution  (e.g.  from  an  embedded  crack  to  a  surface  crack  to  a  through-thickness  crack) 
due  to  intersections  of  growing  cracks  with  surfaces,  as  shown  in  Figure  9,  In  addition, 
some  special  features  of  3D  crack  growth  observed  in  many  real  applications  are  also 
addressed  in  CRACK3D,  such  as  the  disappearance  of  a  crack  front  due  to  local 
separation  of  a  structure,  corner  cracks,  and  cracks  passing  structural  stiffeners.  Figure 
10  shows  the  schematics  of  some  of  the  features  implemented  in  CRACK3D. 

11.10.  Examples 

In  this  section  two  selected  examples  are  used  to  demonstrate  the  reliability  and 
capability  of  CRACK3D  in  the  application  of  simulating  crack  growth  in  three- 
dimensional  components.  In  these  two  examples,  the  effects  of  crack  tunneling  and 
slanting  on  fracture  toughness  are  not  taken  into  consideration.  During  the  process  of 
crack  growth,  the  crack  front  was  assumed  to  be  a  straight  line.  Thus,  the  Mixed-Mode 
CTOD  fracture  criterion  [5]  could  be  employed  in  the  simulations. 

During  the  simulations,  the  direction  of  crack  growth  was  determined  by  the  Mixed- 
Mode  CTOD  fracture  criterion  [5]  and  the  onset  of  crack  growth  was  determined  by  the 
critical  CTOD  of  aluminum  on  the  surface  of  specimens.  When  the  CTOD  at  the 
distance  of  1  mm  behind  the  crack  front  on  the  surface  of  the  specimen  reaches  the 
critical  value  £c=0.08  mm,  then  the  whole  crack  front  extends  0.8mm  in  the  predicted 
direction  of  crack  growth. 

Simulation  of  crack  growth  in  a  plate 

The  geometry  of  the  plate  with  a  single  edge  crack  is  shown  in  Fig.  1 1 .  The  plate  is 
made  of  aluminum  alloy  2024-T3,  and  the  thickness  of  the  plate  is  2  mm.  The  material 
properties  of  aluminum  alloy  2024-T3  are  as  follows.  Young’s  modulus  E  =  71.2 GPa, 
Poisson’s  ratio  v  =  0.3,  initial  yield  stress  oy  =  358  MPa.  The  strain  hardening  curve  of 
2024-T3  is  shown  in  Fig.  12. 

The  plate  underwent  the  Mode  I  displacement  loading,  as  shown  in  Fig.  13.  The  4- 
noded  tetrahedral  elements  were  used  in  the  finite  element  model.  In  the  initial  mesh, 
there  were  5974  elements  and  2155  nodes.  The  deformed  mesh  and  predicted  crack 
path  are  shown  in  Fig.  14.  It  can  be  seen  from  Fig.  14  that,  during  the  process  of  crack 
growth,  the  mesh  in  a  local  region  around  the  crack  front  is  automatically  updated  as  the 
crack  propagates,  and  the  element  density  of  the  mesh  in  all  the  previous  locations  of 
crack  fronts  is  also  updated  during  remeshing  so  that  the  total  number  of  elements  in 
the  finite  element  model  can  be  controlled  to  vary  only  in  a  small  range.  Fig.  15  shows 
the  variation  of  the  number  of  elements  versus  crack  extension  during  the  process  of 
crack  propagation. 

Simulation  of  crack  growth  in  Arcan  specimen 

The  Arcan  test  specimen  and  the  loading  fixture  (composed  of  a  pair  grips)  are 
shown  in  Fig.  16.  The  fixture  is  made  of  15-5PH  stainless  steel  and  has  a  thickness  of 
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12.6  mm.  The  specimen  is  made  of  aluminum  alloy  2024-T3  and  has  a  thickness  of  2.3 
mm.  An  edge  crack  is  introduced  in  the  mid-section  of  the  specimen.  The  initial  crack 
front  is  fabricated  by  the  fatigue  loading.  The  specimen  is  then  attached  to  the  fixture  by 
three  hardened-steel  pins  at  each  end.  The  material  properties  of  15-5PH  stainless 
steel  are  as  follows.  Young’s  modulus  E  =  207 GPa,  Poisson’s  ratio  v  =  0.3,  initial  yield 
stress  Gy  =  M22MPa. 

By  changing  the  angle  of  loading  direction  in  the  fixture,  different  local  mode-mixity 
can  be  obtained.  In  this  section,  only  the  case  of  15  degree  loading  is  provided,  as 
shown  in  Fig.  17.  Considering  that  the  fixture  and  pins  are  relatively  rigid  compared  to 
the  specimen,  the  connection  between  the  fixture  and  specimen  in  the  finite  element 
model  can  be  assumed  to  be  a  rigid  and  continuous  joint.  In  this  simulation  (see  Fig. 
17),  10-noded  tetrahedral  elements  are  employed,  and  there  are  4571  elements  and 
8688  nodes  in  the  initial  mesh.  Also  the  Mixed-Mode  CTOD  fracture  criterion  [5]  is 
employed  in  the  simulation  to  predict  the  onset  and  direction  of  crack  propagation. 

The  profiles  of  deformed  meshes  at  different  loading  steps  are  shown  in  Fig.  16. 
The  comparison  of  predicted  crack  path  to  the  experimentally  measured  crack  path  is 
given  in  Fig.  18.  Fig.  19  shows  the  predicted  load  versus  crack  extension  by 
CRACK3D.  It  is  shown  in  Fig.  19  that  there  is  a  good  agreement  between  the  numerical 
results  and  experimental  results. 

11.11.  Conclusions 

A  finite  element  based,  crack  growth  analysis  and  simulation  code —  CRACK3D 
was  presented.  The  hierarchy  of  the  code  and  some  key  issues  related  to  the  code 
development  and  computational  aspects  of  three-dimensional  crack  growth  simulation 
were  given.  In  the  presented  program  —  CRACK3D,  the  remeshing  technology  is  used 
to  simulate  the  process  of  three-dimensional  crack  growth  in  any  direction  predicted  by 
the  fracture  criteria.  From  the  selected  verification  examples,  it  is  shown  that  CRACK3D 
is  reliable  and  robust.  In  addition  to  being  used  to  simulate  crack  propagation  in  three- 
dimensional  structures  and  components,  it  can  also  be  used  to  develop  and  verify  new 
fracture  criteria  for  brittle  and  ductile  materials. 
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11.12.  Figures 


Figure  1 .  Relation  of  CRACK3D  with  ANSYS  in  the  simulation 


Figure  2.  Hierarchy  of  crack  growth  simulation  code  CRACK3D 
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(d) 


Figure  3.  Schematic  of  crack  front  determination  in  a  plate  with  an  edge  crack 
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Element  quality  y 

(d) 

Figure  4.  3D  mesh  of  a  90°-elbow  pipe  with  a  through-thickness  crack:  (a)  pipe 
geometry  and  crack  location;  (b)  an  overall  view  of  the  mesh;  (c)  a  cross-sectional  view 
along  the  crack  surface;  and  (d)  mesh  quality  distribution. 
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(C)  (d) 


Figure.  5.  3D  mesh  of  a  cube  with  a  semi-elliptical  surface  crack:  (a)  cube  geometry  and 
crack  location;  (b)  an  overall  view  of  the  mesh;  (c)  a  cross-sectional  view  along  the 
crack  surface;  and  (d)  mesh  quality  distribution. 
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Figure.  6.  3D  mesh  of  a  cube  with  an  embedded  elliptical  crack:  (a)  geometry  and  crack 
location;  (b)  an  overall  view  of  the  mesh;  (c)  a  cross-sectional  view  along  the  crack 
surface;  and  (d)  mesh  quality  distribution. 
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(b) 


Figure  7.  A  schematic  oflocal  remeshing  in  CRACK3D,  (a)  mesh  before 
local  remeshing;  (b)  mesh  after  local  remeshing 


Figure  8.  A  new  point  P  in  the  element  Ep  of  the  previous  mesh 
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Figure  9:  Evolution  of  flaw  shapes 
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Figure  10:  Additional  features  during  flaw  evolution 
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Figure  1 1 .  Dimensions  of  a  plate  with  a  single  edge  crack  (Unit:  mm) 
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Figure  12.  Strain  hardening  curves  for  grip  material  and  specimen 
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Figure  15.  Number  of  elements  versus  crack  extension  during  crack 

growth 
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Figure  16.  A  schematic  of  the  Arcan  test  specimen  and  loading  fixture 

(all  dimensions  in  mm) 
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Figure  17.  Initial  finite  element  mesh  and  loading  direction  (0  =15°) 
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(a)  Aa=3.2mm 


(b) 

Aa=1 1.2mm 


(c) 

Aa=19.2mm 


Figure  18.  Crack  propagation  under  mixed-mode  loading 
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Figure  19.  Comparison  of  the  crack  path  (0  =15°) 
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Figure  20.  Experimental  and  numerical  load  versus  crack  extension 


III.  Analysis  and  Simulation  of  Three-Dimensional  Crack  Growth  in  Ductile 
Materials:  Effect  of  Stress  Constraint  on  Crack  Tunneling 

Abstract:  This  work  describes  the  results  of  an  effort  to  understand  crack  tunneling 

events  in  ductile  materials,  which  is  made  possible  by  the  finite  element  based  crack 
growth  simulation  code,  CRACK3D,  where  Section  II  describes  the  underlying 
computational  schemes  developed  for  general  3D  crack  growth  simulation  and  their 
implementation  in  CRACK3D. 

Crack  tunneling  is  a  crack  growth  feature  often  seen  in  stable  tearing  crack 
growth  tests  on  specimens  made  of  ductile  materials  and  containing  through-thickness 
cracks  with  initially  straight  crack  fronts.  As  a  specimen  is  loaded  monotonically,  the 
mid-section  of  the  crack  front  will  advance  first,  which  will  be  followed  by  crack  growth 
along  the  rest  of  the  crack  front,  leading  to  the  formation  of  a  thumb-nail  shaped  crack- 
front  profile.  From  the  viewpoint  of  fracture  mechanics,  crack  tunneling  will  occur  if  the 
operating  fracture  criterion  is  met  first  in  the  mid-section  of  the  crack  front,  which  may 
be  due  to  a  higher  fracture  driving  force  and/or  a  lower  fracture  toughness  in  the  mid¬ 
section.  A  proper  understanding  of  this  fracture  behavior  is  important  to  the 
development  of  a  three-dimensional  fracture  criterion  for  general  stable  tearing  crack 
growth  in  ductile  materials.  In  this  paper,  the  phenomenon  of  crack  tunneling  during 
stable  tearing  crack  growth  in  a  single-edge  crack  specimen  is  investigated  by 
considering  the  effect  of  stress  constraint  on  the  fracture  toughness.  Crack  growth  in  the 
specimen  under  nominally  Mode  I  loading  conditions  is  considered.  In  this  case,  crack 
tunneling  occurs  while  the  initially  flat  crack  surface  (which  is  normal  to  the  specimen’s 
lateral  surfaces)  evolves  into  a  final  slanted  fracture  surface.  A  mixed-mode  CTOD 
fracture  criterion  and  a  custom  three-dimensional  fracture  simulation  code,  CRACK3D, 
are  used  to  analyze  the  tunneling  and  slanting  process  in  the  specimen.  Results  of  this 
investigation  suggest  that  the  critical  CTOD  value  (which  is  the  fracture  toughness)  has 
a  clear  dependence  on  the  crack  front  stress  constraint  (also  called  the  stress  triaxiality, 
which  is  the  ratio  of  the  mean  stress  to  the  von  Mises  effective  stress).  For  simplicity, 
this  dependence  can  be  fitted  to  a  straight  line  within  the  range  of  stress  constraint 
values  found,  with  the  toughness  decreasing  as  the  constraint  increases.  It  is  found  that 
crack  tunneling  in  this  case  is  mainly  the  result  of  a  higher  stress  constraint  (hence  a 
lower  fracture  toughness)  in  the  mid-section  of  the  crack  front.  Details  of  the  crack 
growth  simulation  and  other  findings  of  this  study  will  also  be  presented. 

111.1.  Introduction 

Prediction  of  ductile  fracture  in  metallic  materials  has  been  an  important  subject 
of  fracture  mechanics  research.  A  key  focus  of  this  research  effort  has  been  the 
development  of  fracture  criteria  for  determining  the  onset  and/or  direction  of  crack 
growth.  Currently,  several  approaches  have  been  proposed  to  characterize  the  process 
of  stable  crack  growth  (e.g.  J-Integral  [1],  J-A2  [2],  CTOD  [3,  4]  or  CTOA  [5]  and  so  on). 

Experimental  results  (e.g.  [6-8])  show  that  crack  tunneling  is  a  common  crack 
growth  feature  in  stable  tearing  fracture  in  specimens  made  of  ductile  materials.  For  a 
specimen  containing  a  through-thickness  crack  with  an  initially  straight  crack  front,  as 
the  a  specimen  is  loaded  monotonically,  the  mid-section  of  the  crack  front  will  advance 
first  and  then  the  rest  of  the  crack  front  will  grow,  leading  to  the  formation  of  thumb-nail 
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shaped  crack-front  profiles.  In  addition  to  tunneling,  slant  crack  growth  may  also  occur 
in  ductile  materials.  Stable  tearing  tests  performed  on  specimens  made  of  AL2024-T3 
sheets  indicate  [6,  7]  that  slanting  occurs  when  the  specimens  are  in  the  LT  orientation 
whereas  it  does  not  occur  when  the  specimens  are  in  the  TL  orientation.  On  the  other 
hand,  specimens  made  of  AL  2024-T351  display  both  flat  and  slant  crack  growth  in  LT 
orientation  [8].  However,  what  is  consistent  from  these  tests  is  that  crack  tunneling 
always  occurs  regardless  of  specimen  orientation.  However,  tunneling  in  flat  fracture  is 
usually  more  severe  than  that  in  slant  fracture. 

From  the  viewpoint  of  fracture  mechanics,  the  phenomenon  of  crack  tunneling 
may  be  explained  using  a  fracture  criterion  if  a  non-uniform  fracture  driving  force  and/or 
a  non-uniform  fracture  toughness  exists  along  the  crack  front.  Specifically,  tunneling  will 
occur  if  the  operating  fracture  criterion  is  met  first  in  the  mid-section  of  the  crack  front, 
either  due  to  a  higher  fracture  driving  force  and/or  lower  fracture  toughness  in  the  mid¬ 
section.  A  proper  understanding  of  this  fracture  behavior  is  important  to  the 
development  of  a  three-dimensional  fracture  criterion  for  general  stable  tearing  crack 
growth  in  ductile  materials. 

The  present  work  investigates  crack  tunneling  during  stable  tearing  crack  growth 
in  a  single-edge  crack  specimen  [6,  7],  with  a  focus  on  the  understanding  of  the  effect  of 
stress  constraint  on  fracture  toughness  [9,  10],  based  on  a  CTOD  fracture  criterion  [3, 
4].  It  is  organized  as  follows.  For  reference  purposes,  Section  2  gives  a  brief  description 
of  the  experimental  results  considered  in  this  study.  Section  3  introduces  the  finite 
element  model  and  computational  approach  for  simulating  crack  front  evolution 
observed  in  the  tests.  Section  4  presents  the  distributions  of  CTOD  and  stress 
constraint  along  different  crack  fronts  based  on  the  simulation  results  and  establishes  a 
correlation  between  CTOD  toughness  values  and  stress  constraint  values  in  terms  of  a 
linear  equation.  As  an  application/verification,  this  equation  is  then  used  in  Section  5  as 
part  of  the  CTOD  criterion  to  enable  the  simulation  of  the  stable  tearing  tests  to  predict 
crack  tunneling  profiles  during  crack  growth.  The  predicted  results  are  compared  with 
experimental  measurements.  Section  6  presents  the  conclusions  from  this  study. 

111.2.  Crack  tunneling  measurement 

Figure  1  provides  a  schematic  of  the  single-edge  crack  specimen  geometry  [6]. 
The  specimens  were  machined  from  rolled  2.3mm-thick  sheets  made  of  aluminum  alloy 
2024-T3.  The  specimens  were  then  fatigue  pre-cracked  in  the  LT  orientation  so  that  the 
initial  crack/width  ratio,  a/w=0.0833.  To  assess  the  amount  of  crack  tunneling,  an 
experimental  procedure  described  in  [11]  was  used  to  acquire  crack  front  profile 
measurement  data  on  the  fracture  surface  at  various  levels  of  loading,  as  shown  in 
Figure  2  for  stable  tearing  crack  growth  in  a  specimen  under  remote  Mode  I  loading 
conditions.  The  corresponding  applied  load  for  each  crack  front  profile  is  also  included 
in  the  figure.  It  is  noted  that,  after  the  initiation  of  crack  growth,  the  fracture  surface 
undergoes  a  transition  from  an  initially  flat  crack  plane  to  a  steady  slanted  crack  plane 
as  the  crack  front  advances. 
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III. 3.  Finite  element  model  and  computational  approach 

The  finite  element  method  is  used  to  analyze  the  stable  tearing  crack  growth 
experiments  described  above  in  order  to  understand  the  crack  tunneling  phenomenon. 
To  carry  out  such  analyses,  a  finite  element  code  capable  of  simulating  general  three- 
dimensional  (3D)  crack  growth  in  elastic-plastic  solids  is  required.  In  this  study,  the 
custom  code  CRACK3D  developed  at  the  University  of  South  Carolina  is  used. 
Preprocessing  (for  generating  the  initial  finite  element  mesh)  and  post-processing  (for 
analyzing  stress  and  deformation  state  at  a  particular  stage  of  crack  growth)  are 
performed  using  the  commercial  code  ANSYS,  through  interface  options  in  CRACK3D. 

Two  simulation  options  are  available  in  CRACK3D.  In  the  first  option  (the  nodal 
release  option),  the  crack  front  is  made  to  advance  along  a  prescribed  path,  which  is 
accomplished  through  the  release  of  nodal  pairs  (the  two  nodes  in  a  pair  are  initially  tied 
together  by  rigid  springs)  along  the  crack  path  when  a  certain  condition  (e.g.  when  a 
critical  load  from  a  test  is  reached)  is  satisfied.  This  option  is  useful  for  analyzing 
fracture  tests.  In  the  second  option  (the  local  remeshing  option),  the  crack  front  position 
is  not  prescribed  but  is  predicted  by  a  fracture  criterion  (the  CTOD  criterion  in  this  study) 
and  a  user-specified  region  around  the  new  crack  front  is  remeshed  as  the  crack  front 
grows. 

Figure  3  shows  frontal  planar  views  of  a  3D  finite  element  mesh  used  in 
analyzing  the  stable  tearing  tests  based  on  the  nodal  release  option,  where  Fig.  3(a)  is 
for  the  entire  problem  domain  and  Fig.  (b)  a  local  view  of  the  mesh  around  the  crack 
front  (note  that  the  fatigue  crack  front  has  extended  from  the  initial  notch  into  the  region 
on  the  right).  The  mesh  consists  of  8,917  ten-node  tetrahedral  elements  with  14,315 
nodes.  The  minimal  element  size  around  the  crack  front  is  0.2  mm.  A  remote  Mode  I 
load  was  applied  in  terms  of  a  monotonically  increasing  displacement,  so  that  nodes  at 
the  bottom  edge  of  the  specimen  were  constrained  with  zero  displacements  in  all 
directions,  and  nodes  at  the  specimen’s  top  edge  were  made  to  move  in  the  y-direction 
(vertical  direction)  (while  displacements  in  the  x  and  z  directions  were  held  to  zero). 

Mesh  convergence  analysis  has  been  performed  by  bisecting  each  edge  of 
elements  in  3D  space.  The  study  of  mesh  convergence  here  is  focused  mainly  on  two 
parameters,  the  stress  constraint  and  the  crack  opening  displacement  (COD)  (note  that 
values  of  COD  at  points  not  far  from  the  crack  front  are  called  CTOD  in  this  paper). 
Figures  4  and  5  show  comparisons  of  the  numerical  results  obtained  based  on  finite 
element  mesh  I  (as  shown  in  Fig.  3  with  element  size  of  0.2mm  around  the  crack  front) 
and  finite  element  mesh  II,  which  is  a  refined  mesh  with  element  size  of  0.1mm  around 
the  crack  front.  The  values  of  COD  in  Fig.  4  are  computed  on  the  front  surface  of  the 
specimen,  while  the  values  of  the  stress  constraint  in  Fig.  5  are  computed  on  the  mid¬ 
plane  of  the  specimen.  The  reason  for  choosing  the  mid-plane  of  the  specimen  for  this 
evaluation  is  that  the  gradient  of  the  stress  constraint  ahead  the  crack  front  on  the  mid¬ 
plane  of  the  specimen  is  much  larger  than  that  on  the  free  surfaces  of  the  specimen. 

It  can  be  found  based  on  HHR  fields  [12,  13]  that  both  the  mean  stress  and  the 
von  Mises  effective  stress  have  the  same  singularity  in  terms  of  the  distance  to  crack 
tip.  Therefore,  the  stress  constraint  (which  is  the  ratio  of  mean  stress  and  von  Mises 
effective  stress)  is  expected  to  be  independent  of  the  distance  to  the  crack  tip  in  a  near¬ 
tip  region  under  small  scale  yielding  (SSY)  and  HRR-dominance  conditions.  Indeed 
under  plane  strain  conditions  [14]  the  stress  constraint  is  nearly  constant  along  the 
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radial  direction  in  the  region  around  crack  tip.  However,  neither  the  plane  strain 
condition  nor  the  plane  stress  condition  can  be  exactly  satisfied  in  real-life  cracked 
structures.  Due  to  variations  of  field  singularities  along  three-dimensional  crack  fronts, 
the  stress  constraint  is  expected  to  depend  somewhat  on  the  radial  distance  to  the 
crack  front  even  though  a  strong  near-tip  singularity  is  not  observed  (see  Fig.  5).  To 
alleviate  this  radial  dependence,  an  average  stress  constraint  value  is  taken  in  this 
paper,  which  is  obtained  by  integrating  the  stress  constraint  over  a  small  distance 
ahead  of  the  crack  fronts. 

To  evaluate  fracture  parameter  values  around  the  actual  3D  crack  fronts,  the 
experimentally  measured  fracture  surface  and  crack  front  profiles  associated  with 
different  loading  levels  (see  Fig.  2)  were  built  into  the  finite  element  mesh  to  be  used  in 
the  nodal  release  approach.  It  is  noted  that  the  crack  front  profiles  used  in  the  finite 
element  model  have  been  slightly  smoothed.  Figure  6  shows  the  profiles  of  smoothed 
crack  fronts  in  the  finite  element  model.  Due  to  slant  crack  growth,  the  fracture  surface 
is  not  flat.  As  such,  Fig.  6  is  only  a  projected  section  view  of  the  actually  slant  fracture 
surface  (viewed  from  the  positive  y-direction;  see  Fig.  1).  The  profiles  numbered  1,  2,  4, 
and  6  are  for  crack  fronts  measured  based  on  fatigue  striation  marks  during  interrupted 
crack  growth  tests,  and  the  profiles  numbered  3  and  5  are  interpolated  from  profiles  2 
and  4,  and  4  and  6,  respectively.  Since  the  x-y  plane  coincides  with  the  specimen’s  mid¬ 
plane  (where  z  is  zero),  the  back  surface  of  the  specimen  corresponds  to  negative  z 
coordinates  and  the  front  surface  corresponds  to  positive  z  coordinates.  A  typical  3D 
crack  front  profile  is  shown  in  Fig.  7. 

It  is  noted  that  crack  front  #1  is  the  initial  fatigue  crack  front  and  has  slight 
tunneling  in  the  middle.  As  the  specimen  is  loaded  monotonically  to  the  critical  load  for 
the  initiation  of  crack  growth,  the  middle  region  of  the  specimen  grows  first  while  the 
regions  in  the  specimen’s  front  and  back  surfaces  remain  stationary.  As  the  load  is 
further  increased,  the  crack  front  advances  from  crack  front  #1  to  crack  front  #2.  The 
flaw  eventually  advances  to  crack  front  #6,  with  an  increasing  amount  of  crack 
tunneling.  The  maximum  load  required  for  continued  crack  growth  occurred  at  crack 
front  #6,  after  which  the  required  load  began  to  decrease  gradually. 

CRACK3D  is  used  to  analyze  the  crack  tunneling  and  slanting  process  in  the 
specimen.  With  respect  to  the  fracture  criterion  used  during  the  simulation,  the 
measured  critical  load  corresponding  to  each  crack  front  is  used  to  control  the  crack 
growth  at  different  crack  length.  The  total  CTOD  (defined  as  the  vector  magnitude  of  its 
opening,  shearing  and  tearing  components)  and  stress  constraint  (defined  as  the  ratio 
of  the  mean  stress  to  the  von  Mises  effective  stress)  with  regard  to  the  normal  plane  at 
each  node  on  crack  front  are  evaluated  at  the  critical  instant  just  before  crack 
propagation. 

In  a  CTOD  based  fracture  criterion  [3,  4],  the  driving  force  is  the  total  CTOD, 
which  is  measured  at  a  fixed  distance  behind  the  crack  front.  In  this  study,  CTOD  was 
strictly  computed  at  a  distance  of  0.5  mm  behind  the  crack  front  (along  a  line  normal  to 
the  crack  front).  Since  it  is  impractical  to  use  an  extremely  refined  mesh  and  an 
extremely  small  distance  behind  the  crack  front  for  computing  CTOD,  and  since  severe 
crack  tunneling  creates  a  length  scale  that  may  not  be  sufficient  large  compared  to  0.5 
mm,  the  computed  CTOD  value  is  expected  to  be  less  accurate  in  the  middle  region  of 
a  crack  front  than  away  from  that  region.  The  reason  is  that,  due  to  severe  crack 
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tunneling  in  the  middle  region,  a  point  at  0.5  mm  (or  another  small  distance)  behind  the 
crack  front  along  a  line  normal  to  the  crack  front  may  be  too  close  to  other  parts  of  the 
crack  front,  making  this  point  inappropriate  for  calculating  CTOD.  To  alleviate  this 
problem,  it  is  chosen  in  the  calculation  to  compute  a  nominal  CTOD  value,  when  the 
situation  noted  above  occurs,  by  using  the  crack-front  CTOA  (crack-tip  opening  angle) 
value  to  extrapolate  to  a  distance  of  0.5  mm  behind  the  crack  front. 

It  is  important  to  note  that,  in  this  analysis  phase  of  the  study,  crack  growth 
simulations  were  performed  using  the  nodal  release  option,  with  the  experimentally 
measured  crack  front  profile  positions  and  the  corresponding  loads.  The  CTOD  criterion 
was  not  applied  even  though  CTOD  variations  along  the  crack  fronts  were  evaluated. 

111.4.  Effect  of  stress  constraint  on  critical  CTOD 

The  main  results  of  interests  from  the  finite  element  simulations  described  in  the 
preceding  section  are  the  variations  of  CTOD  and  stress  constraints  along  the 
measured  crack  fronts,  which  will  be  used  to  establish  a  correlation  between  CTOD 
fracture  toughness  and  stress  constraint.  Since  the  CTOD  values  along  a  crack  front 
correspond  to  the  critical  load  that  causes  growth  of  the  crack  front,  they  are  critical 
CTOD  values  and  equal  to  the  corresponding  CTOD  fracture  toughness  values  along 
the  crack  front. 

Figures  8,  9,  10  provide  the  variations  of  the  total  (combined)  CTOD  and  CTOD 
components  for  Mode  I  (opening),  Mode  II  (shearing)  and  Mode  III  (tearing)  along  crack 
fronts  #4,  #5,  and  #6,  respectively,  at  a  distance  of  0.5  mm  behind  the  crack  front  (along 
a  line  normal  to  the  crack  front).  In  the  figures,  the  through-thickness  value  refers  to  the 
z  coordinate  value  along  the  crack  front.  Results  for  crack  fronts  in  the  early  stage  of 
crack  growth  are  not  shown  because  these  crack  fronts  have  growth  only  in  the  middle 
section  and  do  not  provide  reliable  critical  CTOD  values  (note  that  CTOD  values  along 
parts  of  the  crack  front  that  are  not  at  the  impending  moment  of  growth  are  not  critical 
values  and  do  not  equal  to  CTOD  fracture  toughness  values  there). 

Two  important  observations  can  be  made  from  these  figures.  First,  due  to  crack 
tunneling  and  slanting,  a  perfect  symmetry  about  the  specimen’s  mid-plane,  a  feature 
expected  in  Model  I  crack  growth,  is  seen  to  disappear  as  crack  grows,  which  lead  to  a 
truly  three-dimensional  mixed-mode  CTOD  distribution  along  the  crack  front,  especially 
in  the  middle  region  of  the  crack  fronts.  The  crack  front  is  Mode  I  dominant  only  in  the 
middle  region.  Second,  the  total  CTOD  value  (which  represents  the  CTOD-based 
fracture  toughness  during  crack  growth)  is  not  a  constant  along  the  crack  front — it  is 
lower  in  the  middle  region  than  near  the  specimen  surfaces.  As  discussed  later,  this 
variation  in  fact  reflects  the  dependence  of  CTOD  toughness  on  the  stress  constraint 
Am- 

Alternatively,  variations  of  the  total  CTOD  along  the  crack  fronts  can  be  plotted  in 
one  figure,  as  shown  in  Fig.  11.  It  is  seen  that  the  trend  shown  in  the  variations  along 
different  crack  fronts  is  basically  the  same.  This  observation  will  be  utilized 
subsequently  to  relate  CTOD  variation  with  stress  constraint  variation  along  the  crack 
fronts.  Before  this  is  done,  it  must  be  pointed  out  that  CTOD  variations  shown  so  far  are 
strictly  computed  at  a  distance  of  0.5  mm  behind  the  crack  front.  Because  of  severe 
tunneling  (as  discussed  earlier),  the  computed  CTOD  values  are  less  accurate  in  the 
middle  region  of  a  crack  front  than  away  from  that  region.  It  is  suggested  earlier  that 
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improved  CTOD  values  in  the  middle  region  can  be  obtained  by  using  the  crack-front 
CTOA  value  to  extrapolate  to  a  distance  of  0.5  mm  behind  the  crack  front.  This  seems 
to  be  a  viable  approach  for  crack  fronts  after  a  certain  amount  of  crack  growth,  so  that 
crack  blunting  usually  seen  at  the  early  stage  of  crack  growth  can  be  avoided.  Based  on 
this  approach,  the  total  CTOD  variations  in  Fig.  11  are  updated  and  are  shown  in  Fig. 
12. 

Shown  in  Fig.  13  are  the  variations  of  Am  along  crack  fronts  #3,  #4,  #5,  and  #6. 
The  constraint  value  at  each  crack  front  point  is  computed  based  on  an  integrated 
average  from  the  crack  front  to  0.3  mm  ahead  of  the  crack  front,  in  the  direction  normal 
to  the  crack  front,  within  the  plane  of  the  crack  surface.  It  must  be  pointed  out  that,  at 
and  very  near  the  specimen’s  front  and  back  surfaces,  a  distance  of  0.3mm  ahead  of  a 
crack  front  may  go  outside  the  specimen  domain.  In  this  case,  we  use  the  constraint 
value  of  the  next  interior  crack  front  point  to  approximate  the  constraint  value  for  the 
current  crack  front  point.  Similar  to  the  total  CTOD  variations  in  Fig.  12,  a  common  trend 
in  the  stress  constraint  variations  can  also  be  seen. 

To  reduce  scatter  in  the  subsequent  data  reduction  for  possible  correlation 
between  CTOD  toughness  and  stress  constraint  Am,  the  variations  in  Figs.  12  and  13 
need  to  be  smoothed.  To  this  end,  an  averaged  total  CTOD  variation  is  obtained  from 
Fig.  12  and  an  averaged  stress  constraint  variation  is  obtained  from  Fig.  13.  Now,  a 
correlation  between  CTOD  and  Am  can  be  obtained  by  plotting  all  (CTOD,  Am)  pairs 
from  the  same  crack  front  locations,  and  this  correlation  is  shown  in  Fig.  14.  It  is  seen 
that  the  CTOD-based  fracture  toughness  value  decreases  as  the  value  of  the  stress 
constraint  increases.  For  simplicity,  a  straight  line  can  be  used  to  fit  the  CTOD-Am  data 
points  within  the  range  of  values  for  Am,  which  can  be  written  as 

CTOD=0. 0932-0. 031 2  Am  (1) 

Since  the  CTOD  and  Am  values  are  extracted  from  simulation  results  of  actual 
fracture  tests,  with  experimentally  recorded  crack  fronts  and  load  values  during  stable 
crack  growth,  there  is  strong  reason  to  believe  that  Eq.  (1)  may  represent  the 
dependence  of  the  CTOD  based  fracture  toughness  on  stress  constraint  Am  for  the 
material  in  concern,  aluminum  alloy  2024-T3.  Eq.  (1)  suggests  that  the  CTOD  fracture 
toughness  value  is  lower  when  the  stress  constrain  is  higher. 

If  validated  by  further  studies,  the  simple  linear  fit  described  by  Eq.  (1)  or  a  more 
accurate  curve  fit  can  be  treated  as  a  material  property  and  can  be  extremely  useful  in 
three-dimensional  fracture  mechanics  applications  where  constraint  effects  on  fracture 
toughness  are  important. 

III. 5.  Prediction  of  crack  tunneling 

The  relationship  between  the  CTOD  toughness  and  stress  constrain  can  be  used 
to  predict  crack  tunneling  based  on  the  CTOD  fracture  criterion.  To  demonstrate  this 
application,  the  remote  Mode  I  stable  tearing  fracture  test  described  in  previous  sections 
is  now  simulated.  Two  types  of  simulations  using  CRACK3D  are  performed. 

In  the  first  case,  the  nodal  release  option  is  used  with  a  coarse  mesh  (to  save 
computation  time)  with  design  similar  to  that  in  Fig.  3.  The  mesh  consists  of  3,308  ten- 
node  tetrahedral  elements  with  6,169  nodes.  The  minimal  element  size  around  the 
crack  front  is  0.4mm.  The  difference  between  the  simulation  here  and  that  in  section  4  is 
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that,  there  are  no  specified  crack  front  positions  on  the  crack  surface  except  for  the 
initial  fatigue  crack  front  (see  Fig.  15).  The  crack  front  shape  is  predicted  by  the  CTOD 
criterion  and  Eq.  (1 )  even  though  the  crack  front  is  constrained  to  grow  on  the  measured 
slant  fracture  surface.  Since  the  nodal  release  option  is  used,  the  elements  are  fixed.  As 
a  result,  some  mesh  dependence  in  the  predicted  crack  front  positions  is  expected  (see 
Fig.  16).  The  dependence  can  be  reduced  if  finer  meshes  are  used  (a  coarse  mesh  is 
used  here).  Alternatively,  the  crack  front  profile  can  be  smoothed  based  on  predicted 
crack  growth  amounts  at  difference  points  along  the  crack  front. 

In  the  second  case,  instead  of  using  a  fixed  mesh  and  the  nodal  release  option, 
the  local  remeshing  option  is  used  so  that  the  mesh  in  a  user-specified  region  around 
the  current  crack  front  is  remeshed  each  time  the  crack  front  is  predicted  to  grow. 
Again,  the  CTOD  criterion  and  the  relationship  between  the  CTOD  toughness  and 
stress  constraint  are  employed  to  predict  the  crack  front  profile  although  the  crack  front 
is  still  constrained  to  grow  on  the  measured  slant  fracture  surface. 

All  simulations  are  conducted  using  CRACK3D.  In  the  simulations,  the  CTOD 
criterion  is  evaluated  node  by  node  along  the  current  crack  front.  A  crack-front  point  will 
advance  along  the  measured  slant  fracture  surface  when  the  CTOD  value  at  that  point 
reaches  the  critical  value  defined  by  the  stress  constraint  Am  ahead  of  the  crack  front 
according  to  Eq.  (1). 

The  predicted  crack  front  profiles  (without  smoothing)  during  stable  crack  growth 
are  shown  in  Fig.  16  (with  nodal  release)  and  Fig.  17  (with  local  remeshing).  In  order  to 
make  the  different  crack  fronts  displayed  clearly,  some  of  crack  fronts  in  (Fig.  16)  are 
displayed  in  dotted  lines.  From  both  figures  it  was  found  that  the  crack  starts  to  grow  at 
the  mid-section  of  specimen  first,  and  crack  tunneling  happens  as  the  crack  propagates. 
The  shape  of  crack  front  is  consistent  with  the  measured  crack  fronts.  The  depth  of 
crack  tunnel  on  the  crack  path  increases  in  the  early  stage  of  crack  growth,  and 
decreases  after  the  crack  extension  on  the  surface  of  specimen  is  about  1mm,  which  is 
consistent  with  the  measured  results. 


Comparisons  of  predicted  and  measured  crack  tunneling  variations  during  crack 
growth  can  be  made  by  the  use  of  a  non-dimensional  crack  tunneling  depth  parameter, 
as  defined  below  and  in  Fig.  18  (a).  Let  Aai  be  the  amount  of  crack  extension  on  the 
specimen’s  back  surface,  Aa2  the  amount  of  crack  extension  on  the  specimen’s  front 
surface,  and  Aac  the  amount  of  crack  extension  on  the  mid-plane  of  the  specimen.  The 
average  crack  extension  on  the  specimen’s  front  and  back  surfaces  is 


A  a,  +  A  a? 

A  a*  =  — 5 - £ 

s  2 


(2) 


Then  the  non-dimensional  crack  front  tunneling  depth,  T,  is  defined  as 


j._  A  ac  -  A  as 
B 


(3) 


where  B  is  the  thickness  of  the  specimen.  Figure  18  (b)  shows  the  variations  of  T  based 
on  crack-front  profile  measurements  and  simulation  predictions  using  the  nodal  release 
and  local  re-meshing  options  of  CRACK3D.  A  good  agreement  is  observed  in  both 
cases. 
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111. 6.  Conclusions 

In  this  study  the  crack  tunneling  phenomenon  is  investigated  by  considering 
stable  tearing  crack  growth  in  a  single-edge  crack  specimen  made  of  AL  2024-T3. 
Although  the  specimen  was  loaded  remotely  under  Mode  I  conditions,  severe  crack 
tunneling  and  slant  fracture  were  both  present  during  crack  growth.  This  test  was 
analyzed  using  a  custom  three-dimensional  crack  growth  simulation  code,  CRACK3D. 

Results  of  this  study  show  that  the  CTOD-based  fracture  toughness  for  the 
ductile  material  considered  is  a  function  of  stress  constraint  at  the  crack  front.  Based  on 
the  experimentally  recorded  crack  front  profiles  and  the  corresponding  critical  loads,  the 
simulation  results  suggest  that  the  relationship  between  the  CTOD  toughness  and  the 
stress  constraint,  for  AL  2024-T3  and  within  the  stress  constraint  values,  can  be  fitted  to 
a  simple  linear  curve,  such  that  a  higher  stress  constraint  value  would  lead  to  a  lower 
CTOD  toughness  value. 

Based  on  this  study,  crack  tunneling  in  the  specimen  may  be  interpreted  as  the 
result  of  a  lower  CTOD  toughness  (due  to  higher  stress  constraint)  at  the  mid-section  of 
the  crack  front  and  a  higher  CTOD  toughness  (due  to  lower  stress  constraint)  near  the 
specimen’s  front  and  back  surfaces.  Good  comparisons  of  the  predicted  and  measured 
crack  front  profiles  seem  to  confirm  the  result  of  a  previous  study  [10]  that  stress 
constraint  is  a  key  parameter  in  ductile  failure  and  fracture  criteria. 
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III. 8.  Figures 


Fig.  1  Schematic  of  a  single-edge  crack  test  specimen  (all  dimensions  in  mm). 


Fig.  2.  Crack  front  profiles  from  a  Mode  I  stable  tearing  crack  growth  test. 
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(a) 


(b) 


.  3  Frontal  views  of  the  finite  element  mesh:  (a)  mesh  for  the  entire  problem  domain 
and  (b)  zoomed-in  mesh  around  the  crack  front. 
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Front 


Fig.  6  Crack  front  profiles  used  in  stable  tearing  crack  growth  analyses 


Fig.  7  A  3D  crack  front  profile  from  a  stable  tearing  crack  growth  test. 
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Fig.  8  Variation  of  CTOD  and  its  components  along  crack  front  #4 


Through  thickness  (mm) 

Fig.  9  Variation  of  CTOD  and  its  components  along  crack  front  #5 
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Fig.  10  Variation  of  CTOD  and  its  components  along  crack  front  #6 
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Fig.  1 1  Variation  of  total  CTOD  along  crack  fronts  #3,  #4,  #5,  #6  at  0.5  mm 
behind  the  crack  front. 
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Total  CTOD  (mm) 


Through  thickness  (mm) 

Fig.  12  Variation  of  total  CTOD  along  crack  fronts  #3,  #4,  #5,  #6  based  on  improved 
CTOD  calculations  in  the  middle  region  of  the  crack  front. 


Fig.  13  Variation  of  stress  constraint  Am  along  crack  fronts  #3,  #4,  #5,  #6 
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Fig.  14  Correlation  between  the  critical  CTOD  and  stress  constraint  A* 


Initial  crack 


Fig.  15  Finite  element  mesh  on  the  fracture  surface 
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Fig.  17  Comparisons  of  predicted  (solid  lines,  with  local  remeshing)  and  measured  (symbols) 
crack-front  profiles. 
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Fig.  18  Crack  front  tunneling  during  crack  growth:  (a)  definition  of  a  non-dimensional  crack 
tunneling  depth,  (b)  comparisons  of  measured  and  predicted  crack  tunneling  depth 
variations  with  crack  growth. 
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IV.  A  combined  experimental  and  finite  element  study  of  crack  slanting  in  a 
ductile  material  under  mixed-mode  loading 

Abstract.  Slant  fracture  is  widely  observed  during  crack  growth  in  thin  sheet  specimens 
made  of  ductile  materials,  providing  a  good  case  for  investigating  three-dimensional 
criteria  for  mixed-mode  ductile  fracture.  To  gain  an  understanding  of  slant  fracture 
events  and  to  provide  insight  for  establishing  a  slant  fracture  criterion,  stable  tearing 
fracture  experiments  on  combined  tension-torsion  (nominal  mixed-mode  I/Ill)  specimens 
and  nominal  Mode  I  Arcan  specimens  made  of  Al  2024-T3  are  analyzed  using  the  finite 
element  method  under  three-dimensional  conditions.  Two  types  of  finite  element  models 
are  considered  for  the  study  of  slant  fracture:  (a)  combined  tension-torsion  specimens 
containing  stationary,  flat  and  slant  cracks  subject  to  loads  corresponding  to  the  onset 
of  crack  growth,  and  (b)  stable  tearing  crack  growth  with  slanting  in  a  nominal  Mode  I 
Arcan  specimen.  Analysis  results  reveal  that  there  exists  a  strong  correlation  between 
the  direction  of  the  maximum  effective  plastic  strain  ahead  of  a  crack  front  and  the 
orientation  direction  of  slant  fracture.  In  particular,  it  is  observed  that  (a)  at  the  onset  of 
crack  growth,  the  angular  position  of  the  maximum  effective  plastic  strain  around  the 
crack  front  serves  as  a  good  indicator  for  the  slant  fracture  surface  orientation  during 
subsequent  crack  growth;  and  (b)  during  stable  tearing  crack  growth  with  a  flat-to-slant 
transition,  the  crack  growth  path  on  each  section  plane  through  the  thickness  of  a 
specimen  coincides  with  the  angular  position  of  the  maximum  effective  plastic  strain 
around  the  crack  front.  The  results  of  this  study  suggest  that  the  effective  plastic  strain 
may  be  used  as  a  fracture  parameter  for  predicting  the  fracture  surface  path  of  stable 
crack  growth  with  slant  fracture. 

IV.1.  Introduction 

A  phenomenon  often  observed  in  stable  tearing  crack  growth  experiments  on  thin 
plate  specimens  made  of  ductile  materials  (e.g.  aluminum  alloys)  is  that  an  initially  flat 
crack  tends  to  turn  and  grow  as  a  slant  crack  after  a  short  flat-to-slant  transition.  Figure 
la  shows  a  typical  crack  surface  involving  the  development  of  slant  fracture  under 
nominal  Mode  I  loading  conditions.  This  behavior  has  been  reported  widely  in  fatigue 
experiments  (Rickerby  and  Fenici,  1984;  Zuidema  and  Blaauw,  1988;  Zehnder,  2000), 
stable  tearing  tests  (Meyn  et  al.,  1989;  Narasimhan  et  al.,  1992;  Amstutz,  1995; 
Mahmoud  and  Lease,  2003),  and  dynamic  crack  growth  tests  (Krafft  et  al.,  1961). 

Experiments  on  compact  tension  or  C(T)  specimens  with  initial,  slant  cracks  have 
been  carried  out  for  different  ductile  materials,  for  instance,  steels  (Kumar  and  Hirth, 
1991)  and  aluminum  alloys  (Manoharan,  1997).  Different  slant  cracks  were  used  to 
obtain  different  Mode  III  and  Mode  I  ratios  under  remote  Mode  I  loading  conditions. 
These  studies  indicate  that  slant  cracks  tend  to  lower  the  critical  value  of  the  total 
mixed-mode  J-integral,  Jmc.  To  this  end,  it  is  worth  noting  that  experiments  on  middle- 
crack-tension  or  M(T)  specimens  made  of  Al  2024-T351,  conducted  by  Dawicke  and 
Sutton  (1994),  showed  that  the  critical  crack-tip  opening  angle  (CTOA)  value  measured 
on  a  specimen  surface  is  higher  in  the  case  of  a  flat  crack  than  in  the  case  of  a  slant 
crack. 
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Crack-front  shapes  in  Mode  I  loaded  Al  2024-T351  plate  specimens,  recorded  at 
different  stages  of  crack  growth,  have  been  compared  by  James  and  Newman  (2003). 
The  optically  measured  fracture  crack-front  shapes  revealed  that  specimens 
experiencing  a  flat-to-slant  transition  have  much  lower  tunneling  magnitudes  than  those 
with  only  flat  fracture  surfaces.  Here  the  tunneling  parameter  defined  by  Dawicke  and 
Sutton  (1994),  which  is  the  normalized  difference  between  the  maximum  interior  crack 
extension  and  the  crack  extension  on  the  specimen  surface,  is  used  to  define  tunneling 
magnitude. 

Under  pure  tension  and  combined  tension-torsion  (mixed-mode  I/Ill)  loading 
conditions,  experiments  conducted  by  Sutton  et  al.  (2001)  on  single-edge  crack 
specimens  made  of  aluminum  alloy  2024-T3  suggest  that  the  slant  angle,  which  is  the 
angle  between  the  slant  plane  and  the  original  flat  crack  plane,  has  a  one-to-one 
correspondence  with  the  ratio  of  a  nominal  shear  stress,  Sr,  due  to  a  remote  torque 
load  T,  to  a  nominal  tensile  stress,  SP,  due  to  a  remote  tension  load  P  (Fig.  1b).  For 
each  fixed  loading  ratio,  the  slant  angle  initially  increases  with  the  amount  of  crack 
extension  during  a  flat-to-slant  transition  and  eventually  achieves,  more  or  less,  a 
steady-state  value.  It  is  found  that  the  slant  angle  decreases  as  the  amount  of  torsion 
loading  increases.  When  only  tension  is  applied,  the  slant  angle  is  -38°,  which  is  the 
largest  slant  angle.  It  is  noted  that  the  same  slant  angle  is  also  observed  in  Arcan 
experiments,  also  on  Al  2024-T3  specimens,  by  Amstutz  (1995)  under  remote  Mode  I 
loading  conditions. 

Analyses  of  fracture  surfaces  of  M(T)  specimens  made  of  Al  2024  have  been 
performed  by  Bron  et  al.  (2004).  Comparisons  of  flat  and  slant  fracture  surfaces 
revealed  that  the  formation  of  the  slant  plane  is  related  to  shear  band  localization.  An 
earlier  study  by  Randolph  and  Piascik  (1995)  also  showed  that  shear  lips  were  formed 
along  the  specimen’s  free  surfaces  at  the  beginning  of  the  stable  tearing  process  and 
then  they  grew  to  form  a  slant  crack. 

Based  on  continuous  damage  mechanics,  a  numerical  study  of  the  flat-to-slant 
transition  was  carried  out  by  Besson  et  al.  (2001)  using  the  Rousselier  model 
(Rousselier,  1987).  It  was  assumed  that  crack  growth  was  controlled  by  porosity  and 
that  a  material  failed  when  porosity  reached  a  critical  value.  The  authors  made  a 
qualitative  prediction  of  the  flat  to  slant  transition  by  changing  certain  model  parameters 
such  as  the  initial  plastic  strain  for  void  nucleation.  However,  this  model  overestimated 
the  structural  response  and  underestimated  the  crack  growth  rate.  Gullerud  and  his 
coauthors  (Gullerud  et  al.,  1999)  argued  that  the  lack  of  agreement  with  experimental 
evidence  was  due  to  the  fact  that  this  model  required  a  high  constraint  level  to  drive  the 
damage  process.  One  approach  to  obtain  an  improved  prediction  (Steglich  and  Brocks, 
1998)  is  to  use  small  finite  elements  with  element  size  on  the  order  of  the  inter-particle 
spacing,  which  is  about  10-20  pm  for  Al  2024,  a  very  small  number. 

Analytical  and  numerical  studies  of  slant  fracture  using  traditional  fracture 
mechanics  concepts  have  been  very  limited  in  the  literature.  Under  linearly  elastic 
conditions,  the  crack  tip  stress  fields  for  plates  loaded  in  tension  and  out-of-plane  shear 
was  derived  by  Zehnder  et  al.  (2000)  based  on  superposition  of  the  stresses  from  in¬ 
plane  and  out-of-plane  loads.  Through  an  analysis  of  the  crack-tip  stress  fields,  the 
authors  attempted  to  explain  why  the  crack  took  the  slant  orientation. 
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To  shed  light  on  the  underlying  mechanics  of  slant  fracture,  three-dimensional 
finite  element  analyses  of  selected  Al  2024-T3  specimen  geometries  with  initial  flat  and 
slant  cracks  under  nominal  Mode  I  loading  condition  were  carried  out  by  Maghoub  et  al. 
(2003).  For  the  specimen  with  a  slant  crack,  the  slant  angle  was  set  to  38°,  which  is  the 
value  from  a  nominal  Mode  I  experiment  of  the  same  material  (Amstutz,  1995). 
Comparisons  of  deformation  and  stress  fields  around  the  flat  and  slant  crack  fronts 
revealed  that,  among  other  observations,  the  von  Mises  effective  stress  ahead  of  a  slant 
crack  is  more  uniform  through  the  plate  thickness  and  is  higher  than  observed  in  the 
case  of  a  flat  crack  undergoing  the  same  loading  conditions.  The  findings  in  Mahgoub  et 
al.  (2003)  are  consistent  with  the  observations  in  Bron  et  al.,  (2004).  The  authors 
indicated  that  the  effective  stress  and  other  stress  field  features  promote  shear  fracture, 
which  is  believed  to  be  related  to  the  formation  of  slant  fracture. 

Slant  fracture  usually  involves  mixed-mode  I/I  1 1  loading  conditions.  Due  to  the 
complexity  of  general  mixed-mode  I/Ill  problems,  analytical  solutions  of  the  3D  crack 
front  stress  and  deformation  fields  do  not  exist.  However,  under  more  idealized 
conditions,  such  as  linear  elasticity,  small  scale  yielding,  Mode  I  or  Mode  III  dominance, 
and/or  with  negligible  field  variations  in  the  plate  thickness  direction,  a  limited  number  of 
studies  have  been  published,  such  as  those  by  Pan  (1990),  Pan  and  Shih  (1990),  Hui 
and  Zehnder  (1993)  and  Gao  and  Shih  (1998). 

Currently  there  are  no  fracture  criteria  for  quantitative  prediction  of  slant  fracture  and 
many  slant  fracture  issues  remain  unresolved.  For  example,  under  combined  tension- 
torsion  loads  (Sutton  et  al.,  2001),  one  may  ask  whether  is  it  possible  (a)  to  predict  the 
experimentally  observed  orientation  of  the  slant  fracture  surface  that  tends  to  create 
contact  and  interference  between  the  two  crack  surfaces  instead  of  opening  and 
clearance  between  the  surfaces,  and  (b)  to  predict  the  experimentally  observed 
dependence  of  the  stable  slant  angle  on  the  torsion/tension  loading  ratio? 

The  purpose  of  the  current  work  is  to  gain  a  further  understanding  of  the  slant 
fracture  phenomenon  and  to  provide  insight  for  establishing  a  slant  fracture  criterion  for 
quantitative  prediction  of  stable  tearing  crack  growth  with  slant  fracture.  To  this  end, 
three-dimensional  (3D)  elastic-plastic  finite  element  analyses  of  stable  tearing 
experiments  involving  slant  fracture  are  carried  out.  The  focus  will  be  on  two  types  of 
experiments:  (a)  combined  tension-torsion  (nominal  mixed-mode  I/Ill)  experiment  and 
(b)  a  nominal  Mode  I  Arcan  experiment.  The  specimens  are  made  of  Al  2024-T3 
material.  It  will  be  shown  that  the  angular  position  of  the  maximum  effective  plastic 
strain  around  the  crack  front,  which  coincides  with  the  direction  of  the  maximum  extent 
of  the  effective  plastic  strain  contours  around  the  crack  front,  correlates  strongly  with  the 
orientation  of  the  slant  fracture  surface  and  with  the  crack  growth  paths  on  section 
planes  through  the  thickness  of  a  specimen.  Based  on  the  analysis  results,  it  will  be 
argued  that  the  effective  plastic  strain  may  be  used  as  a  fracture  parameter  for 
predicting  the  fracture  surface  orientation  during  slant  fracture.  It  is  believed  that  the 
findings  of  this  study  can  serve  as  a  basis  for  formulating  a  fracture  criterion  for 
predicting  mixed-mode  ductile  fracture  events  under  3D  conditions. 

Subsequent  sections  are  arranged  as  follows.  In  Section  2,  3D  finite  element 
models  of  the  combined  tension-torsion  experiments  (Sutton  et  al.,  2001)  and  the  Mode 
I  Arcan  (Amstutz,  1995)  on  Al  2024-T3  specimens  are  described.  In  Section  3,  crack- 
front  strain  fields  for  the  combined  tension-torsion  experiments  at  the  onset  of  crack 
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growth,  obtained  under  3D,  elastic-plastic  and  large  deformation,  are  presented  and 
their  implications  investigated.  Section  4  provides  the  results  of  a  3D  stable  tearing 
crack  growth  simulation  along  an  experimentally  observed  3D  fracture  surface,  of  a 
Mode  I  Arcan  experiment  with  crack  slanting.  Finally,  a  summary  of  the  main  findings  of 
this  study  and  concluding  remarks  are  given  in  Section  5. 

IV.2.  3D  Finite  Element  Models 

IV. 2.1  Combined  Tension-Torsion  Experiments 

Combined  tension-torsion  experiments  on  Al  2024-T3  specimens  have  been 
carried  out  by  Sutton  et  al.  (2001)  to  investigate  stable  tearing  crack  growth  under 
remote  mixed-mode  I/Ill  loading  conditions.  These  experiments  are  analyzed  in  this 
study  to  gain  an  understanding  of  possible  correlations  between  the  crack-front 
deformation  states  at  the  onset  of  crack  growth  and  the  slant  fracture  surface  orientation 
during  subsequent  crack  growth. 

Figure  2  shows  the  in-plane  geometry  of  the  plate  specimen  in  a  global  Cartesian 
coordinate  system  (the  specimen  thickness  is  2.3  mm).  The  specimen  was  fatigue  pre¬ 
cracked  in  the  LT  orientation  (hence  the  tension  loading  is  in  the  rolling  or  longitudinal 
direction  and  the  crack  is  parallel  to  the  transverse  direction),  with  an  initial  a/w=0.083, 
where  a  is  the  pre-crack  length  and  w  the  width  of  the  plate.  The  pin  holes  near  the  top 
and  bottom  specimen  surfaces  were  used  to  clamp  the  specimen  to  a  tension-torsion 
loading  fixture  consisting  of  a  pair  of  loading  tangs,  grip  plates,  and  a  backing  plate  for 
alignment  and  installation  (which  was  detached  prior  to  the  start  of  an  experiment). 

Three  loading  cases,  ST/SP=0.00,  ST/SP=1.66  and  ST/SP=6.64,  are  considered  in 
this  study,  where  St  is  a  nominal  value  of  the  maximum  shear  stress  due  to  the  applied 
torque  T  and  SP  is  a  nominal  value  of  the  average  normal  stress  due  to  the  tensile  load 
P,  and  they  are  given  by 

j  ST  =3T  (t 2  w)'1 
I  Sp  =  P  (t  w)"1 

where  t  is  the  specimen  thickness  and  w  is  the  specimen  width.  It  is  noted  that 
ST/SP=0.00  corresponds  to  a  Mode  I  loading  condition. 

For  illustration  purposes  in  this  and  subsequent  sections,  Figure  3  shows  a 
schematic  drawing  of  a  plate  specimen  with  a  slant  crack  and  a  crack-front  rectangular 
Cartesian  coordinate  system  in  which  the  y  axis  is  along  the  tensile  loading  direction 
and  the  x-y  plane  coincides  with  the  specimen’s  mid-plane.  Three  of  the  specimen’s 
through-thickness  section  planes,  namely  the  front  surface,  the  mid-plane,  and  the  back 
surface,  which  correspond  to,  respectively,  z=0.5f,  0,  and  -0.5f,  are  also  noted  in  the 
figure.  When  the  slant  angle  is  zero  the  specimen  contains  an  initially  flat  crack; 
otherwise  the  specimen  contains  an  initially  slant  crack.  A  crack-front  local  polar 
coordinate  system  (r,  9)  can  be  associated  with  a  through-thickness  section  plane  (e.g. 
the  front  surface,  the  mid-plane,  and  the  back  surface),  as  shown  in  Fig.  3  for  the  mid¬ 
plane,  with  r= 0  at  the  crack  front  and  0=0°  parallel  to  the  positive  x  direction. 

Actual  material  properties  for  the  specimen  and  the  loading  fixture  are  employed. 
The  aluminum  specimen  has  a  Young’s  modulus  of  71.7  GPa  and  an  initial  yield  stress 
of  344.8  MPa,  while  the  loading  fixture,  which  is  made  of  steel  15-5PH  (a  high-strength 
stainless  steel),  has  a  Young’s  modulus  of  207  GPa  and  an  initial  yield  stress  of  1 ,724 
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MPa.  A  Poisson’s  ratio  of  0.3  is  taken  for  both  the  steel  and  the  aluminum  alloy.  Both 
materials  exhibit  strain-hardening  behavior  and  the  dependence  of  the  materials’  flow 
stress  on  the  effective  plastic  strain  in  uniaxial  tension  is  shown  in  Fig.  4.  The  materials 
are  assumed  to  obey  the  J2  flow  theory  of  plasticity. 

A  finite  element  representation  of  the  tension-torsion  specimen  with  a  flat  crack, 
with  consideration  of  the  four  grip  plates  in  the  fixture-specimen  connection  region,  is 
illustrated  in  Fig.  5.  The  mesh  shown  in  Fig  5  is  a  converged  mesh  used  for  analysis.  In 
particular,  Fig.  5a  provides  an  overall  view  of  the  finite  element  model.  Figure  5b  is  a 
local  in-plane  view  of  the  focused  mesh  around  the  crack  front;  the  crack  is  not  visible 
but  it  lies  to  the  left  of  the  center  point  of  the  mesh.  In  the  radial  direction,  the  element 
size  decreases  towards  the  crack  front.  In  the  circular  direction,  forty  rings  of  elements 
are  distributed  uniformly  around  the  crack  front.  Figure  5c  shows  a  three-dimensional 
view  of  the  mesh  in  a  near-crack-front  region  with  a  radius  of  4.5  mm.  In  the  thickness 
direction,  there  are  sixteen  layers  of  elements,  with  a  decreasing  layer  thickness  from 
the  plate  mid-plane  towards  the  front  and  back  surfaces  (the  element  size  on  the 
surface  is  one  fifth  of  that  in  the  mid-plane).  This  mesh  consists  of  26,304  twenty-node 
brick  elements  with  113,985  nodes  (25,072  elements  in  the  specimen  region  and  1,232 
elements  in  the  fixture  region).  Because  of  large  out-of-plane  deformation  in  the 
combined  tension-torsion  experiments,  it  is  necessary  to  apply  a  large  deformation  finite 
element  formulation.  The  general-purpose  finite  element  code  ANSYS  was  employed 
with  the  updated  Langrangian  method. 

The  specimen  is  loaded  as  follows.  The  bottom  surface  of  the  lower 
fixture/specimen  region  is  held  fixed  while  the  top  surface  of  the  upper  fixture/specimen 
region  is  loaded  according  to  the  conditions  specified  in  Table  1,  which  correspond  to 
the  critical  moment  of  the  onset  of  stable  tearing  fracture.  Specifically,  the  upper  grip  is 
rotated  by  an  experimentally  measured  torque  T  (which  is  applied  through  a  set  of 
statically  equivalent  point  forces)  and  displaced  vertically  by  Uy  that  results  in  a  reaction 
Fy  in  agreement  with  the  measured  value.  The  use  of  a  large  deformation  formulation 
and  the  applied  loading  process  described  above  have  been  shown  (Lan  et  al.,  2004)  to 
be  effective  for  modeling  the  thin  tension-torsion  specimen.  It  is  noted  that  this  particular 
loading  approach  for  applying  accurate  boundary  conditions  is  required  since  reliable 
experimental  data  was  available  only  for  the  tensile  reaction  force,  Fy,  the  applied 
torque,  T  (about  the  y  axis)  and  the  resulting  torsion  angle  </>  (rotation  about  the  y  axis). 
The  good  agreement  shown  in  Table  1  between  the  predicted  and  measured  specimen 
response  in  terms  of  Fy  and  <f>  demonstrates  that  the  large  deformation  option  and 
applied  finite  element  loading  conditions  are  appropriate. 

IV.  2. 2  Mode  I A  ream  Experiment 

As  a  special  case  of  the  mixed-mode  I/ll  Arcan  stable  tearing  crack  growth  experiments 
performed  by  Amstutz  (1995)  on  Al  2024-T3  specimens,  the  Mode  I  experiment  on  a 
specimen  with  the  LT  orientation  (in  which  the  initial  crack  orientation  is  perpendicular  to 
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Table  1  Finite  element  loading  conditions,  predicted  specimen  response,  and  test 
data  for  three  tension-torsion  test  cases 


Test  Case 

Top 

Surface 

Conditions 

Predicted  Response 

Test  Data 

ST/Sp=0.0 

0 

Uy=0.48 

mm 

Reaction  Fy=37.56  KN 
Rotation  angle  <f>  =  0° 

Reaction  Fy=37.19  KN 
Rotation  angle  </>  =  0° 

Sj/Sp 

=1.66 

Uy=0.28 

mm 

T=32.06 

N*m 

Reaction  Fy=24.85  KN 
Rotation  angle  <j>  =  10.96° 

Reaction  Fy=25.24  KN 
Rotation  angle  (j>  =  11.83° 

St/Sp 

=6.64 

Uy=- 

0.03mm 

T=42.91 

N*m 

Reaction  Fy=8.91  KN 
Rotation  angle  </)  =  17.24° 

Reaction  Fy=8.43  KN 
Rotation  angle  (/>  =  17.83° 

the  material’s  rolling  direction)  provides  a  good  case  for  investigating  the  slant  fracture 
phenomenon  and  is  chosen  for  detailed  analysis  in  the  current  work.  In  this  Mode  I 
Arcan  experiment,  as  the  specimen  containing  an  initially  flat  fatigue  pre-crack  is  loaded 
under  remote  Mode  I  loading  conditions,  the  flat  crack  first  experiences  a  flat-to-slant 
transition  growth  period  and  then  grows  as  a  slant  crack  with  an  almost  constant  slant 
angle. 

Figure  6  shows  the  in-plane  geometry  of  the  Arcan  fixture  and  specimen,  in 
which  the  specimen  has  a  thickness  of  2.3  mm  and  the  fixture  (made  of  15-5PH 
stainless  steel)  has  a  thickness  of  19  mm.  The  locations  of  the  pinholes  on  the  outer 
edge  of  the  fixture  provide  a  range  of  values  for  the  loading  angle <f>  where 

^  =  0°  corresponds  to  Mode-1  loading  and  ^  =  90°  corresponds  to  Mode-ll  loading.  The 
three  interior  holes  on  each  of  the  fixture  halves  are  used  to  attach  the  specimen  to  the 
fixture. 

The  crack  growth  event  in  the  Arcan  Mode  I  specimen  is  modeled  using  a  nodal 
release  technique  in  which  the  prescribed  crack  growth  path  is  based  on  the  measured 
fracture  surface.  Shown  in  Fig.  7  is  a  converged  3D  finite  element  mesh  of  the 
fixture/specimen  system  used  for  the  crack  growth  simulation.  Since  the  fixture  and 
connecting  pins  in  the  Arcan  experiments  are  relatively  rigid  compared  to  the  aluminum 
specimen,  the  fixture-specimen  system  is  treated  as  one  continuous  solid  with  two 
regions  of  different  thicknesses  and  material  properties  (Deng  and  Newman,  1999). 
This  mesh  consists  of  12,101  ten-node  tetrahedral  elements  with  21,040  nodes.  The 
smallest  element  size  along  the  crack  surface  is  *  0.2  mm. 

The  stable  tearing  crack  growth  process  (see  Fig.  7c),  starting  from  a  flat  fatigue 
pre-crack,  following  a  flat-slant  transition  region,  and  settling  into  a  slant  crack  growth 
region  with  a  slant  angle  *38°,  is  simulated  using  CRACK3D,  which  is  a  finite  element 
code  developed  at  the  University  of  South  Carolina  for  3D  crack  growth  simulation  using 
nodal  release  and  local  re-meshing  options  (see  Zuo  et  al.,  2004,  2005). 

Crack  growth  along  the  measured  fracture  surface  in  Fig.  7c  is  achieved  through 
the  advancement  of  an  assumed  straight  crack  front  (thus  crack  tunneling  is  not 
modeled)  via  a  nodal  release  procedure,  in  which  the  crack  front  nodes  (originally  tied 
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together  by  rigid  springs)  are  separated  (released).  Nodal  release  is  performed  when, 
based  on  the  mixed-mode  crack  tip  opening  displacement  (CTOD)  criterion  (Ma  et  al., 
1999;  Sutton  et  al.,  2000),  the  generalized  CTOD  value  on  the  front  surface  of  the 
specimen  and  at  a  fixed  distance  behind  the  crack  tip  reaches  an  experimentally 
determined  critical  value.  In  this  study,  an  experimentally  determined  critical  CTOD 
value  of  0.095  mm  at  1 .0  mm  behind  the  crack  tip  (which  corresponds  to  a  critical  CTOA 
value  of  5.44°)  is  used.  Unlike  the  combined  tension-torsion  experiments,  large  out-of¬ 
plane  deformation  is  not  present  in  Arcan  specimens,  and  as  such,  it  is  found  that  only 
the  small-deformation  option  in  CRACK3D  is  needed  to  properly  simulate  the  crack 
growth  process. 

IV.3.  Results  for  the  combined  tension-torsion  experiments 

Tension-torsion  specimen  models  with  stationary  flat  and  stationary  slant  cracks 
have  been  analyzed  using  loading  conditions  corresponding  to  the  onset  of  crack 
growth.  For  models  containing  a  flat  crack,  the  finite  element  analysis  results  will  show 
that  there  is  a  strong  correlation  between  the  direction  of  the  maximum  effective  plastic 
strain  at  the  onset  of  crack  growth  and  the  subsequent  slanting  direction.  This 
correlation  will  be  further  demonstrated  by  results  from  models  with  a  slant  crack,  which 
suggests  that  the  effective  plastic  strain  may  serve  as  a  good  indicator  for  predicting 
crack  growth  with  slant  fracture. 

IV.3. 1  Tension-Torsion  Models  with  a  Flat  Crack 

Effective  plastic  strain  contours 

The  effective  plastic  strain  contours  on  the  specimen’s  front  surface  (z=t/2)  just  before 
the  onset  of  crack  growth  is  shown  in  Fig.  8  for  three  tension-torsion  loading  cases  (with 
torsion/tension  ratios  ST/SP=0.00,  1 .66  and  6.64).  Besides  the  observation  that  the  size 
of  the  contour  extent  increases  with  the  torsion/tension  ratio,  the  results  also  reveal  that 
the  orientation  angle  (relative  to  the  x-axis)  of  the  maximum  extent  of  the  effective 
plastic  strain  contours  on  the  front  surface  decreases  with  the  torsion/tension  ratio. 
Specifically,  for  the  contour  level  of  0.0001,  the  maximum  contour  extent  for  the  pure 
tension  case  (St/Sp=0.00)  occurs  at  0  *  +/-46°  relative  to  the  x-axis,  while  those  for  the 
mixed  tension-torsion  cases  the  maximum  0  «  37°  (Sj/SP=1.66)  and  0  *  36° 
(St/Sp=6.64)  on  the  front  surface  occurs.  On  the  back  surface,  the  angles  for 
ST/SP=1 .66  and  6.64  are  the  same  as  those  on  the  front  surface  but  with  negative  signs, 
which  is  due  to  an  anti-symmetry  between  the  contour  orientations  on  the  front  and 
back  surfaces  (compare,  for  example,  Fig.  8b  and  Fig.  9  for  the  case  of  St/Sp=1  .66).  On 
the  specimen’s  mid-plane,  the  effective  plastic  strain  contours  are  symmetric  about  the 
x-axis  and  the  maximum  extent  occurs  more  or  less  along  the  x-axis  direction  (see  Fig. 
10). 

To  make  connections  with  experimental  observations  (Sutton  et  al.,  2001),  it  is 
noted  that  the  orientation  of  the  slant  fracture  surface  in  a  tension-torsion  experiment  is 
such  that  the  presence  of  a  torsion  load  component  will  tend  to  create  crack  surface 
contact  and  interference  instead  of  opening  and  clearance.  Specifically,  for  a  mixed 
tension-torsion  loading  case,  the  subsequent  slant  orientation  after  the  onset  of  crack 
growth  is  such  that  the  crack  path  on  the  model  specimen’s  front  surface  will  tilt 
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upwards  while  the  crack  path  on  the  back  surface  will  go  downwards.  This  is  entirely 
consistent  with  the  orientation  of  the  maximum  extent  of  the  effective  plastic  strain 
contours,  as  shown  in  Fig.  8b  (front  surface)  and  Fig.  9  (back  surface)  for  the  case  of 
St/Sp=1.66.  Furthermore,  experimental  results  (Fig.  1b)  show  that  the  slant  angle 
decreases  as  the  torsion/tension  load  ratio  increases  (the  slant  angle  is  0  *  38°,  33.5° 
and  31°,  respectively,  for  ST/Sp=0.00,  1.66  and  6.64).  This  is  also  consistent  with  the 
numerical  results  that  the  orientation  angle  of  the  maximum  extent  of  the  effective 
plastic  strain  contours  decreases  on  the  front  surface  with  the  torsion/tension  ratio. 

Angular  variations  of  the  effective  plastic  strain  field 

In  order  to  reinforce  the  observation  that  there  is  a  strong  correlation  between  the 
effective  plastic  strain  distribution  and  the  orientation  direction  of  slant  fracture,  the 
angular  variations  of  the  effective  plastic  strain  around  the  crack  front  is  examined. 
These  angular  variations  are  given  for  a  through-thickness  section  plane  (e.g.  the  front 
surface,  the  mid-plane,  and  the  back  surface)  of  the  specimen,  in  terms  of  the  local 
polar  coordinate  0  (see  Fig.  3),  with  r=  0  at  the  crack  front  and  0=0°  parallel  to  the 
positive  x  direction. 

Figure  11  shows  the  angular  variations  of  the  effective  plastic  strain  ep 
(normalized  by  the  initial  yield  strain  £o)  along  a  circular  path  of  radius  r=1.50mm,  on  the 
mid-plane  and  the  front  and  back  surfaces,  for  three  tension-torsion  loading  cases.  It  is 
seen  from  Fig.  11a  that  the  effective  plastic  strain  on  the  front  surface  reaches  a 
maximum  at  0  *  +/- 45°,  34°  and  30°,  respectively,  for  ST/SP=0.00,  1.66  and  6.64.  Hence 
the  absolute  value  of  the  angle  at  which  the  maximum  effective  plastic  strain  occurs  on 
the  front  surface  decreases  as  torsion  increases.  This  trend  is  strongly  indicative  of  the 
trend  observed  experimentally  for  the  slant  angle,  namely,  that  the  absolute  value  of  the 
slant  angle  decreases  as  torsion  increases.  Specifically,  experimental  results  showed 
that,  after  a  short  flat-to-slant  transition,  the  slant  angle  in  tension-torsion  experiments  is 
0  «  38°,  33.5°  and  31°,  respectively  for  Sj/SP=0.00,  ST/SP=1 .66  and  St/Sp=6.64.  For  the 
case  of  ST/Sp=0.00,  the  slant  angle  may  be  0  *  38°  or  -38°,  but  for  cases  with  a  torsion 
loading  component  (e.g.  ST/SP=1.66  and  ST/Sp=6.64),  the  slant  angle  is  always  positive, 
which  confirms  the  trend  observed  for  the  predicted  maximum  effective  plastic  strain 
angles  on  the  front  surface. 

Consistent  with  the  positions  of  the  maximum  effective  plastic  strain  on  the  front 
surface,  those  on  the  back  surface  have  the  same  but  opposite  angular  values  (see  Fig. 
11c).  As  such,  if  these  angles  are  indicators  of  where  the  crack  path  will  tend  to  turn  on 
the  front  and  back  surfaces,  these  indicators  will  correctly  predict  the  qualitative  trends 
exhibited  by  the  subsequent  slant  fracture  orientations  seen  from  the  experimental 
results:  (a)  when  torsion  is  present,  the  slant  orientation  will  tend  to  created  crack 
surface  contact  and  interference  instead  of  crack  opening  and  clearance;  and  (b)  the 
slant  angle  decreases  as  the  torsion/tension  ratio  increases.  The  observation  (see  Fig. 
11b)  that  the  angle  at  which  the  maximum  effective  plastic  strain  occurs  on  the 
specimen’s  mid-plane  is  0°  (for  Sj/Sp=0.00  and  Sj/Sp=1.66)  or  somewhere  around  0° 
(for  ST/SP=6.64)  is  also  consistent  with  the  predicted  overall  slant  orientation  through 
the  specimen  thickness. 
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The  above  comparisons  demonstrate  again  that  there  is  a  strong  correlation  between 
the  effective  plastic  strain  distribution  and  the  direction  of  crack  growth  with  slant 
fracture.  This  effective  plastic  strain  argument  will  be  discussed  again  in  Section  3.2. 

Other  stress  and  strain  field  variations 

For  completeness  and  comparison,  other  relevant  field  variations  at  the  onset  of 
crack  growth  are  included  here.  First,  the  radial  variations  of  the  stress  constraint  Om/oe 
(which  is  the  ratio  between  the  mean  stress  om  and  the  von  Mises  effective  stress  cre) 
and  the  effective  plastic  strain  sp  will  be  presented,  and  then  the  angular  variations  of 
om,  oe  and  Om/ae  will  be  given. 

Figure  12  shows  the  radial  variations  of  the  constraint  Om/oe  along  the  direction 
0=0°  ahead  of  a  flat  crack  on  the  front  surface,  mid-plane,  and  the  back  surfaces,  with 
three  tension-torsion  loading  cases.  It  is  seen  that,  in  all  loading  cases,  the  mid-plane 
constraint  value  (Fig.  12b)  rises  quickly  as  the  crack  front  is  approached  and,  near  the 
crack  front,  it  is  much  higher  than  those  on  the  front  and  back  surfaces  (Figs.  12  a,  c, 
which  show  almost  identical  variations,  suggesting  some  kind  of  constraint  symmetry 
about  the  mid-plane).  It  is  also  observed  that  the  constraint  value  decreases  as  the 
torsion/tension  ratio  increases. 

Figure  13  describes  the  radial  variations  of  the  effective  plastic  strain  (normalized 
by  the  initial  yield  strain  e0)  along  0=0°  on  the  front  surface  (the  front  and  back  surfaces 
have  identical  variations)  and  on  the  mid-plane.  It  is  seen  that,  close  to  the  crack  front 
(when  r<0.3mm),  £p  increases  rapidly  as  the  crack  front  is  approached  and  as  the 
torsion/tension  ratio  increases,  which  is  especially  true  on  the  front  and  back  surfaces. 
Figures  14,  15,  and  16  gives  the  angular  variations  of  the  von  Mises  effective  stress,  the 
mean  stress,  and  the  constraint,  respectively,  along  r=1.50mm  for  the  mid-plane  and  the 
front  and  back  surfaces  for  three  loading  cases.  Based  on  the  data  in  the  figures,  two 
observations  are  noted.  First,  the  effective  stress  (Fig.  14)  is  more  or  less  flat  ahead  of 
the  crack  front,  which  suggests  that  it  may  not  serve  well  as  a  parameter  for  predicting 
the  crack  path.  Second,  the  mean  stress  (Fig.  15)  and  the  constraint  (Fig.  16)  have  a 
peak  value  somewhere  around  0=0.  More  specifically,  for  ST/SP=0.00,  the  peak  value 
occurs  exactly  at  0=0°  on  the  mid-plane  and  front  and  back  surfaces,  and  for 
ST/SP=1.66  and  6.46,  it  occurs  at  0=0°  only  on  the  mid-plane. 

IV. 3. 2  Tension-Torsion  Models  with  a  Slant  Crack 

In  the  preceding  section  it  has  been  argued  that  the  effective  plastic  strain  distribution  at 
the  onset  of  crack  growth  in  tension-torsion  specimens  containing  a  flat  crack  provides  a 
good  indicator  for  the  orientation  of  the  subsequent  slant  fracture  surface.  In  this 
section,  further  evidence  will  be  provided  for  this  argument.  To  this  end,  one 
modification  is  made  to  the  tension-torsion  specimens  in  the  3D  finite  element  analyses 
performed  earlier:  the  flat  crack  is  replaced  with  a  slant  crack  while  all  other  input  data, 
such  as  material  properties,  geometrical  dimensions,  and  boundary  conditions  are  kept 
the  same.  In  doing  so,  it  is  hoped  that  the  following  issues  can  be  clarified:  (1)  If  the 
slant  angle  of  the  stationary  slant  crack  equals  the  experimentally  measured  slant 
angle,  will  the  effective  plastic  strain  distribution  indicate  a  constant  slant  angle  during 
subsequent  crack  growth?  (2)  If  the  slant  angle  is  smaller  than  the  measured  value,  will 
the  effective  plastic  strain  distribution  indicate  an  increase  in  the  slant  angle  during 
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subsequent  crack  growth?  (3)  If  the  slant  angle  is  larger  than  the  measured  value,  will 
the  effective  plastic  strain  distribution  indicate  a  decrease  in  the  slant  angle  during 
subsequent  crack  growth? 

To  address  these  questions,  3D  finite  element  analyses  for  tension-torsion 
specimen  models  containing  stationary  slant  cracks  have  been  carried  out.  The  slant 
crack  models  have  the  same  material,  geometry  and  applied  remote  loading  conditions 
as  those  with  flat  cracks  except  that  the  crack  surface  in  the  slant  crack  model  forms  a 
slant  angle,  a ,  with  that  of  the  flat  crack  surface  (see  Fig.  3).  For  each  of  the  three 
tension-torsion  loading  cases  ST/Sp=0.00,  1.66,  and  6.64  (for  which  the  measured  slant 
angles  are,  respectively,  0  *  38°,  33.5°  and  31°),  five  slant  angles  have  been 
considered:  0°,  10°  and  20°  (which  are  smaller  than  the  measured  slant  angle),  55° 
(which  is  larger  than  the  measured  slant  angle),  and  the  measured  slant  angle  itself. 
The  results  are  presented  below.  Since  conclusions  from  the  results  for  the  three 
tension-torsion  cases  are  the  same,  only  those  for  the  case  of  ST/Sp=1.66  are  reported 
here. 

Angular  variations  of  the  effective  plastic  strain  field 

Figure  17  shows  the  angular  variations  of  the  effective  plastic  strain  along 
r=1.50mm  on  the  front  surface  (Fig.  17a)  and  on  the  middle  plane  (Fig.  17b),  for  the 
case  of  St/Sp=1.66  and  for  the  five  slant  angle  values  (a=0°,  10°,  20°,  33.5°and55°).  It 
is  clearly  seen  that,  when  the  slant  angle  equals  0°,  10°,  or  20°,  which  are  smaller  than 
the  experimentally  observed  slant  angle  of  33.5°,  the  angle  at  which  the  maximum 
effective  plastic  strain  occurs  on  the  front  surface  (Fig.  17a)  is  positive  (it  is  worth  noting 
that  the  corresponding  angle  on  the  back  surface  the  same  value  but  opposite  sign). 
According  to  the  argument  made  in  the  preceding  section,  this  would  suggest  that  in  the 
subsequent  crack  growth,  the  crack  path  on  the  front  surface  will  curve  upwards  while 
that  on  the  back  surface  will  go  downwards,  thus  resulting  an  increase  in  the  slant 
angle.  Furthermore,  it  is  observed  that  the  angle  for  the  maximum  effective  plastic  strain 
will  decreases  as  the  slant  angle  gets  closer  to  the  measured  value.  Hence,  based  on 
the  effective  plastic  strain  argument,  the  predicted  increase  in  the  slant  angle  for 
subsequent  crack  growth  will  be  smaller  when  the  difference  between  the  stationary 
slant  angle  and  the  measured  angle  is  smaller. 

On  the  other  hand,  when  the  slant  angle  is  55°,  which  is  larger  than  the 
measured  slant  angle  of  33.5°  ,  the  angle  at  which  the  maximum  effective  plastic  strain 
occurs  is  negative  on  the  front  service  (and  positive  on  the  back  surface).  Based  on  the 
effective  plastic  strain  argument,  this  would  suggest  that  the  crack  path  in  subsequent 
crack  growth  will  turn  downwards  on  the  front  surface  and  upwards  on  the  back  surface, 
resulting  in  a  decrease  in  the  slant  angle.  It  is  interesting  to  note  that,  on  the  mid-plane 
(Fig.  17b)  the  angle  at  which  the  maximum  effective  plastic  strain  occurs  is  basically 
zero  regardless  of  the  slant  angle,  which  suggests  that  the  crack  path  on  the  mid-plane 
will  stay  on  the  straight  path,  a  result  that  is  consistent  with  experimental  observations. 
Another  important  observation  is  that,  when  the  slant  angle  equals  the  measured  value 
of  33.5° ,  the  angle  at  which  the  maximum  effective  plastic  strain  occurs  is  always  zero 
on  all  surfaces  (see,  e.g.,  Fig.  17),  regardless  of  the  tension-torsion  ratio.  This  is  further 
confirmed  by  Fig.  18,  which  shows  that  the  maximum  effective  plastic  strain  peak  is 
always  located  at  0=0°  on  all  section  planes  through  the  thickness  of  the  specimen. 
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Based  on  the  effective  plastic  strain  argument,  this  set  of  observations  would  suggest 
that  the  crack  path  on  each  section  plane=through  the  specimen  thickness  will  not 
change  direction.  Thus  the  slant  angle  during  subsequent  crack  growth  will  stay 
constant,  which  is  confirmed  by  the  tension-torsion  experimental  results. 

Table  2  Slant-angle  dependence  of  the  angle  at  which  the  maximum  effective  plastic 
strain  occurs  on  the  front  and  back  surfaces  for  three  tension-torsion  loading  cases  (a* 
denotes  the  experimentally  observed  stable  slant  angle  for  the  three  loading  cases). 


ST/S 

p 

Slant  angle  (  Degree) 

0° 

10° 

20° 

* 

a 

55° 

Fron 

t 

Back 

Fron 

t 

Back 

Front 

Bac 

k 

Front 

Bac 

k 

Front 

Back 

0.00 

±45 

±45 

40.5 

-40.5 

31.5 

31.5 

0 

0 

-31.5 

31.5 

1.66 

34 

-34 

27 

-27 

13.5 

13.5 

0 

0 

-40.5 

40.5 

6.64 

30 

-30 

18 

-18 

9 

-9 

0 

0 

-54 

54 

A  summary  of  the  predicted  slant-angle  dependence  of  the  angle  at  which  the  maximum 
effective  plastic  strain  occurs  on  the  front  and  back  surfaces  for  all  three  tension-torsion 
loading  cases  is  given  in  Table  2.  This  table  and  the  analyses  presented  in  the 
preceding  paragraphs  clearly  demonstrate  that  the  effective  plastic  strain  distribution 
provides  a  good  indicator  for  predicting  the  fracture  surface  paths  during  stable  crack 
growth  with  slant  fracture.  This  argument  will  be  further  supported  by  findings  from  a 
finite  element  simulation  of  stable  crack  growth  along  an  experimentally  observed  slant 
fracture  surface  in  a  nominally  Mode  I  loaded  Arcan  specimen,  which  will  be  described 
in  Section  4. 

Other  stress  and  strain  field  variations 

For  comparison,  the  radial  variations  of  constraint  and  effective  plastic  strain  ahead  of 
the  slant  crack  front  (at  the  onset  of  crack  growth)  are  included  for  the  case  of 
St/Sp=1.66.  Figure  19  shows  the  radial  variations  of  the  constraint  along  9=0°  on  the 
front  surface  and  the  mid-plane,  around  slant  cracks  with  various  slant  angles 
(variations  on  the  back  surface  are  not  shown  because  they  are  basically  the  same  as 
those  on  the  front  surface).  It  is  seen  that  the  constraint  level  on  the  front  surface  (Fig. 
19a)  is  always  lower  than  0.6,  while  on  the  mid-plane  (Fig.  19b)  it  rises  quickly  (with  the 
exception  of  a=55°)  as  the  crack  front  is  approached.  Also  on  the  mid-plane,  constraint 
decreases  as  the  slant  angle  increases.  Except  for  the  case  a=55°,  the  constraint  level 
on  the  mid-plane  is  always  greater  than  that  on  the  front  surface  (Fig.  19a). 

Figure  20  presents  the  radial  variations  of  the  effective  plastic  strain  along  0=0° 
around  slant  cracks  with  various  slant  angle;  variations  on  the  back  surface  are  identical 
to  those  on  the  front  surface  (Fig.  20a).  Variations  on  the  mid-plane  (Fig.  20b)  are 
supplemented  with  a  closer  view  near  the  crack  front  (Fig.  20c).  As  can  be  seen  from 
Figs.  20a  and  c,  the  effective  plastic  strain  near  the  crack  front  increases  when  the  slant 
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angle  increases  from  0°  to  33.5°  (which  is  the  measured  slant  angle  for  ST/SP=1.66), 
and  its  value  decreases  when  the  slant  angle  is  55°. 

IV.4.  Results  for  the  Mode  I  Arcan  Experiment 

Results  presented  in  Section  3  are  from  finite  element  models  containing  stationary,  flat 
and  slant  cracks  subjected  to  loading  conditions  corresponding  to  the  onset  of  stable 
tearing  crack  growth.  Effective  plastic  strain  distributions  around  the  stationary  crack 
fronts  are  correlated  to  the  orientation  of  slant  fracture  surface  in  the  subsequent  crack 
growth  events. 

To  further  establish  the  correlation  between  slant  fracture  and  the  effective  plastic  strain 
and  lay  a  solid  foundation  for  developing  a  criterion  for  predicting  slant  fracture,  this 
section  focuses  on  results  from  finite  element  simulations  of  an  actual  stable  tearing 
crack  growth  event  involving  slant  fracture,  as  reported  by  Amstutz  (1995)  on  an  Al 
2024-T3  specimen  for  a  Mode  I  Arcan  experiments,  in  which  an  initially  flat  fatigue  pre¬ 
crack  turns  into  a  slant  crack  as  the  crack  front  advances.  Figure  7  shows  the  finite 
element  model  for  this  specimen.  The  model  contains  a  prescribed  crack  path  based  on 
measurements  of  the  fracture  surface  profile  after  the  experiment.  The  stable  tearing 
crack  growth  process  is  simulated  using  CRACK3D  (see  Zuo  et  al.  2004)  using  a  nodal 
release  option  so  that  the  evolution  of  the  effective  plastic  strain  field  variations  can  be 
captured  and  examined  in  detail  as  the  crack  undergoes  flat-to-slant  transition. 

Figure  21  provides  nine  snapshots  abstracted  from  simulation  movies  of  the 
fracture  event  as  seen  on  the  front  surface,  mid-plane,  and  the  back  surface,  revealing 
the  evolution  of  the  effective  plastic  strain  contours  on  (a)  the  front  surface,  (b)  the  mid¬ 
plane,  and  (c)  the  back  surface,  with  crack  extension  increments  (1)Aa  =  2.0mm, 
(3)Aa  =  4.0  mm  and  (3)  Aa  =  6.0  mm  .  For  example,  Fig.  21. b2  shows  the  effective  plastic 
strain  contours  on  the  mid-plane  when  the  crack  has  advanced  2  mm.  The  white  line  in 
each  figure  represents  the  actual  crack  path  on  each  section  plane,  which  turns 
downwards  on  the  front  surface,  stays  almost  straight  on  the  mid-plane,  and  moves 
upwards  on  the  back  surface.  These  figures  clearly  demonstrate  that,  as  the  crack 
advances,  the  crack  path  on  each  section  plane  basically  follows  the  direction  of  the 
maximum  extent  of  the  effective  plastics  strain  contours. 

Figure  22  presents  the  angular  variations  of  the  effective  plastic  strain  along 
r=1.5  mm  on  the  front  surface,  mid-plane,  and  the  back  surface  for  various  amounts  of 
crack  extension  from  A  a  =  0.00  mm  to  Aa  =  8.00  mm .  Several  important  observations  can 
be  made  from  these  figures.  First,  each  angular  variation  curve  on  the  mid-plane  (Fig. 
22b)  has  only  one  maximum  and  it  always  occurs  at  about  0=0°,  which  is  consistent 
with  the  observation  that  the  crack  path  on  the  mid-plane  is  always  along  the  direction 
of  0=0°  (i.e.  the  crack  path  stays  straight).  Second,  except  at  the  onset  of  crack  growth 
(i.e.  A  a  =  0.00  mm ),  each  angular  variation  on  the  front  and  back  surface  has  only  one 
maximum,  which  occurs  at  an  angle  that  is  negative  or  approximately  0°  on  the  front 
surface  and  is  positive  or  approximately  0°  on  the  back  surface.  Third,  the  maximum 
value  of  the  effective  plastic  strain  increases  as  A  a  increases  and  the  absolute  value  of 
the  angle  associated  with  the  maximum  value  decreases  gradually  and  eventually 
settles  down  to  0°  as  the  crack  experiences  the  flat  to  slant  transition  (the  flat-to-slant 
transition  region  is  about  4  mm,  as  shown  in  Fig.  7c)  and  stays  on  the  stable  slant 
fracture  surface.  That  is,  as  the  crack  grows  the  direction  of  the  maximum  effective 
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plastic  strain  on  each  section  plane  becomes  more  and  more  aligned  with  the  crack 
path  direction.  Stated  more  explicitly,  the  direction  of  the  maximum  effective  plastic 
strain  on  each  section  plane  always  points  towards  the  correct  crack  growth  direction. 

The  second  and  third  observations  from  Fig.  22  are  expressed  quantitatively  in 
Table  3,  which  lists  the  direction  (the  angle)  of  the  maximum  effective  plastic  strain  on 
the  three  section  planes  for  crack  extension  increments  from  Aa  =  0.00  mm 

to  Aa  =  8.00  mm .  The  double  values  (i.e.  0  =  +/-45°)  for  the  front  and  back  surfaces  at 

the  onset  of  crack  growth  are  expected  because  of  a  mathematical  bifurcation  for  crack 
slanting  for  a  specimen  containing  a  perfectly  flat  crack  under  Mode  I  loading 
conditions.  In  practice,  the  actual  crack  slanting  direction  is  affected  by  imperfections. 
For  example,  the  initial  fatigue  pre-crack  is  probably  not  perfectly  flat,  that  is,  the  crack 
surface  may  have  an  initial  slant  angle,  leading  to  a  unique  direction  for  the  maximum 
effective  plastic  strain  at  the  onset  of  crack  growth.  To  see  this,  consider  the  case  in 
which  the  initial  fatigue  crack  surface  has  an  exaggerated  initial  slant  angle  of  -10°. 
Then  a  unique  maximum  effective  plastic  strain  direction  of  about  -36°  (instead  of  0  «  +/- 
45°)  will  be  observed  for  the  front  surface.  Of  course,  if  a  smaller  initial  slant  angle  is 
used,  the  unique  maximum  effective  plastic  strain  direction  will  be  closer  to  0  *  -45°  on 
the  front  surface. 

Table  3  Direction  of  the  maximum  effective  plastic  strain  on 
the  three  section  planes  for  A  a  =  0.00  mm  to  Aa  =  8.00  mm . 


Crack 

extension 

(mm) 

Front 

surface 

(degree) 

Mid-plane 

(degree) 

Back 

surface 

(degree) 

0.00 

±45 

0 

±45 

1.00 

-31.5 

-1 

40.5 

2.00 

-25 

-5 

17 

3.00 

-7 

-5 

9 

4.00 

-4 

1 

3 

5.00 

1 

2 

3 

6.00 

0 

1 

0 

7.00 

0 

0 

0 

8.00 

0 

0 

0 

A  more  visually  convincing  representation  of  the  above  observations  can  be  seen  from 
Fig.  23,  which  plots  the  variation  of  an  angle  6*  with  crack  extension  on  the  front  surface 
(Fig.  23a),  the  mid-plane  (Fig.  23b)  and  the  back  surface  (Fig.  23c).  The  angle  Q*  is 
relative  to  the  x-axis  and  is  positive  if  counterclockwise.  In  the  case  of  experimental 
data,  the  angle  represents  the  crack  growth  direction  (the  kink  angle).  In  the  case  of 
simulation  data,  two  sets  of  values  are  denoted  by  6*\  (a)  the  angular  position  around 
the  crack  tip  (at  r=1.5  mm  away  from  the  crack  tip)  where  the  maximum  sp  occurs,  and 
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(b)  the  direction  with  the  maximum  radial  extent  of  sp  contours.  This  figure  clearly 
illustrates  the  equivalence  between  the  direction  of  the  maximum  extent  of  the  effective 
plastic  strain  contours  and  the  angle  at  which  the  maximum  effective  plastic  strain 
occurs,  and  the  strong  correlation  of  the  maximum  effective  plastic  strain  angle  with  the 
crack  growth  path  on  each  section  plane  and  hence  with  the  overall  orientation  of  the 
slant  fracture  surface.  Taken  together,  these  observations  and  previous  findings  point  to 
the  conclusion  that  the  effective  plastic  strain  ahead  of  a  crack  front  provides  a  good 
indicator  for  predicting  the  direction  of  crack  growth  path  on  each  section  plane  through 
the  thickness  of  a  specimen,  and  hence,  for  predicting  the  overall  orientation  of  the  slant 
fracture  surface. 

IV.5.  Summary  and  concluding  remarks 

The  phenomenon  of  slant  fracture  in  several  combined  tension-torsion  specimens  and 
in  a  nominal  Mode  I  Arcan  specimen,  all  made  of  Al  2024-T3,  has  been  studied  using 
the  finite  element  method. 

For  the  tension-torsion  specimens,  3D  finite  element  models  containing  stationary,  flat 
and  slant  cracks  subject  to  loads  at  the  onset  of  crack  growth  have  been  analyzed 
under  elastic-plastic  and  large  deformation  conditions.  Results  show  that  the  direction  of 
the  maximum  effective  plastic  strain  at  the  onset  of  crack  growth  strongly  correlates  with 
the  orientation  of  the  slant  fracture  surface  during  subsequent  stable  tearing  crack 
growth. 

For  the  nominal  Mode  I  Arcan  specimen,  a  3D  finite  element  model  containing  a 
flat  fatigue  pre-crack  and  a  measured  slant  fracture  surface  for  the  subsequent  stable 
tearing  crack  growth  path  has  been  used  to  simulate  the  actual  slant  fracture  process  in 
the  specimen.  Simulation  results  confirm  the  findings  from  the  tension-torsion 
experimental  analyses  and  further  demonstrates  that  the  direction  of  the  maximum 
effective  plastic  strain  (or  the  direction  of  the  maximum  extent  of  the  effective  plastic 
strain  contours)  provides  a  good  indicator  for  the  direction  of  the  crack  growth  path  on 
section  planes  through  the  thickness  of  the  specimen,  and  hence  with  the  overall 
orientation  of  the  slant  fracture  surface. 

The  findings  of  this  study  suggest  that  the  direction  of  the  maximum  effective 
plastic  strain  may  be  utilized  as  an  effective  parameter  for  predicting  the  slant  fracture 
orientation.  While  these  findings  contribute  to  the  understanding  of  mechanics  issues  in 
slant  fracture  and  may  provide  a  basis  for  formulating  a  3D  mixed-mode  fracture 
criterion  for  predicting  crack  growth  events  in  ductile  materials,  much  work  still  remains 
to  be  done.  For  example,  in  addition  to  the  effective  plastic  strain,  stress  constraint  is 
also  known  to  be  a  key  parameter  in  controlling  ductile  failure  processes  (e.g. 
McClintock,  1968;  Rice  and  Tracey,  1969;  Hancock  and  Mackenzie,  1976;  Walsh  et  al., 
1989;  Bao  and  Wierzbicki,  2004).  An  important  open  issue  is  the  role  of  the  dependence 
of  the  critical  effective  plastic  strain  on  the  stress  constraint  level  in  the  prediction  of 
ductile  crack  growth  events.  It  is  hoped  that  the  findings  presented  in  this  paper  will  offer 
useful  insights  for  further  work  in  this  research  area. 
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IV.7.  Figures 


Fig.  1  (a)  A  typical  slant  fracture  surface  showing  an  initially  flat  crack  surface  and  a  flat-to-slant 
transition  region,  and  (b)  evolution  of  slant  angle  during  combined  tension-torsion  tests  in 
Al  2024-T3  specimens  (Sutton,  et  al.,  2001 ),  as  a  function  of  the  amount  of  crack 


Fig.  2  Combined  tension-torsion  test  specimen  geometry  (all  dimensions  in  mm) 
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Fig.  3  A  schematic  of  a  plate  specimen  showing  a  slant  crack,  a  crack-front 
coordinate  system,  and  a  local  polar  coordinate  system. 
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Fig.  4  Strain  hardening  curves  for  the  specimen  (Al 
2024-T3)  and  fixture  (15-5PH  stainless  steel) 
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(a)  (c) 


Fig.  5  A  converged  finite  element  mesh  for  the  tension-torsion  specimen: 
(a)  a  global  view,  showing  four  grip  plates,  (b)  a  local  in-plane  view  of  the 
focused  mesh  around  the  crack  front,  (c)  a  three-dimensional  view  of  the 
mesh  in  a  near-crack-front  region. 
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Fig.  6  In-plane  geometry  of  the  Arcan  fixture  and 
specimen:  (a)  fixture;  (b)  specimen. 
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(a) 


(b) 


Fig.  7  Finite  element  mesh  for  an  Arcan  LT  specimen  loaded  in  Mode  I 
(viewed  from  the  front):  (a)  the  global  mesh,  (b)  a  local  view  focusing  on 
the  crack  path  region,  and  (c)  a  local  view  showing  different  regions  and 
their  dimensions  of  the  fracture  surface  along  a  measured  flat-slant  crack 
growth  path. 
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Fig.  9  Effective  plastic  strain  contours  on  the  specimen  back  surface, 
for  the  tension-torsion  loading  case  of  Sr/SP=1 .66. 


Fig.  10  Effective  plastic  strain  contours  for  a  strain  level  of  0.10e-3, 
on  the  specimen  middle  plane,  for  three  tension-torsion  loading 
cases. 


72 


-180  -135  -90  -45  0  45  90  135  180 

6  (Degree) 

(c) 

Fig.  1 1  Angular  variations  of  the  effective  plastic  strain  (normalized  by  the 
initial  yield  strain)  along  r=1.50  mm  around  a  flat  crack:  (a)  on  the  front 
surface,  (b)  on  the  middle-plane,  and  (c)  on  the  back  surface. 
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Fig.  12  Radial  variations  of  the  constraint  along  0=0°  ahead  of  a  flat  crack: 
(a)  on  the  front  surface,  (b)  on  the  middle-plane,  and  (c)  on  the  back  surface 
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Fig.  13  Radial  variation  of  the  effective  plastic  strain  (normalized  by  the  initial 
yield  strain)  along  0=0°  ahead  of  a  flat  crack:  (a)  on  the  front  and  back 
surfaces  and  (b)  on  the  mid-plane. 
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Fig.  14  Angular  variations  of  the  von  Mises  effective  stress  along  r=1 .50  mm 
around  a  flat  crack:  (a)  on  the  front  surface,  (b)  on  the  mid-plane,  and  (c)  on  the 
back  surface. 
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Fig.  15  Angular  variations  of  the  mean  stress  along  r=1.50  mm  around  a  flat 
crack:  (a)  on  the  front  surface,  (b)  on  the  mid-plane,  and  (c)  on  the  back  surface. 
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Fig.  16  Angular  variations  of  the  constraint  along  r=1.50  mm  around  a  flat  crack: 
(a)  on  the  front  surface,  (b)  on  the  mid-plane,  and  (c)  on  the  back  surface. 
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Fig.  17  Angular  variations  of  the  effective  plastic  strain  along  r=1.50  mm  around 
slant  cracks  with  various  slant  angles:  (a)  on  the  front  surface  and  (b)  on  the 
mid-plane. 
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Fig.  18  Angular  variations  of  the  effective  plastic  strain  along  r=1 .50  mm  around  a 
slant  crack  that  equals  the  measured  value  on  several  through-thickness  section 
planes. 
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Fig.  20  Radial  variations  of  the  effective  plastic  strain  along  0=0°  around  slant 
cracks  with  various  slant  angles:  (a)  on  the  front  surface,  (b)  on  the  mid-plane,  and 
(c)  a  local  view  of  variations  on  the  mid-plane. 
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Fig.  21  Evolution  of  the  effective  plastic  strain  contours  on  (a)  the  front  surface,  (b) 
the  mid-plane,  and  (c)  the  back  surface,  with  a  crack  extension  amount  of 

(1 )  A  a  =  2.0  mm  ,  (3)  A  a  =  4.0  mm  and  (3)  A  a  =  6.0  mm  . 
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Fig.  22  Angular  variations  of  the  effective  plastic  strain  along  r=1 .5  mm  on  (a)  the 
front  surface,  (b)  the  middle  plane  and  (c)  the  back  surface,  for  various  amounts  of 
crack  extension. 
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Fig.  23  Variation  of  0*  with  crack  extension  on  (a)  the  front  surface,  (b)  the  mid- 
plane  and  (c)  the  back  surface,  where  G*  is  the  crack  growth  direction  (for  test), 
the  angle  at  which  the  maximum  sp  occurs  (for  simulation),  or  the  direction  with 
the  maximum  extent  of  sp  contours  (for  simulation).  The  angle  6*  is  relative  to  the 
x-axis  and  is  positive  if  counterclockwise. 
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